<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Analyze Steady-State Community COBRA Models</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S3 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S4 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S5 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: pre;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: pre;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S6 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S7 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S14 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Analyze Steady-State Community COBRA Models</span></h1><h2  class = 'S1'><span>Author(s): Siu Hung Joshua Chan, Department of Chemical Engineering, The Pennsylvania State University</span></h2><h2  class = 'S1'><span>Reviewer(s): Almut Heinken, Luxembourg Centre for Systems Biomedicine, University of Luxembourg</span></h2><div  class = 'S2'><span style=' font-style: italic;'></span></div><h2  class = 'S1'><span>INTRODUCTION</span></h2><div  class = 'S2'><span>This tutorial demonstrates the use of SteadyCom to analyze a multi-organism COBRA model (e.g., for a microbial community) at a community steady-state [1]. Compared to the direct extension of flux balance analysis (FBA) which simply treats a community model as a multi-compartment model, SteadyCom explicitly introduces the biomass variables to describe the relationships between biomass, biomass production rate, growth rate and fluxes. SteadyCom also assumes the existence of a time-averaged population steady-state for a stable microbial community which in turn implies a time-averaged constant growth rate across all members. SteadyCom is equivalent to the reformulation of the earlier community flux balance analysis (cFBA) [2] with significant computational advantage. SteadyCom computes the maximum community growth rate by solving the follow optimization problem:</span></div><div  class = 'S3'><span texencoding="\begin{array}{ll}
\max &amp; \ \mu\\ \\
\text{s.t.} &amp; \sum\limits_{j\in \textbf{J}^k}S^k_{ij}V^k_j=0, &amp; \forall i \in \textbf{I}^k, k \in \textbf{K}\\
 &amp; LB^k_jX^k\leq V^k_j\leq UB^k_jX^k, &amp; \forall j \in \textbf{J}^k, k \in \textbf{K} \\
&amp; \sum\limits_{k \in \textbf{K}}V^k_{ex(i)} + u^{com}_i\geq 0, &amp; \forall i \in \textbf{I}^{com} \\
&amp; V^k_{biomass} = X^k\mu, &amp; \forall k \in \textbf{K} \\
&amp; \sum\limits_{k \in \textbf{K}}X^k = 1 \\
&amp; X^k,\quad \mu \geq 0, &amp; \forall k \in \textbf{K} \\
&amp; V^k_j \in \Re, &amp; \forall j \in \textbf{J}^k, k \in \textbf{K} 
\end{array}" style="vertical-align:-124px"><img src="" width="280.5" height="259" /></span></div><div  class = 'S2'><span>where </span><span texencoding="S^k_{ij}" style="vertical-align:-8px"><img src="" width="17.5" height="22" /></span><span> is the stoichiometry of metabolite </span><span style=' font-style: italic;'>i</span><span> in reaction </span><span style=' font-style: italic;'>j</span><span> for organism </span><span style=' font-style: italic;'>k</span><span>, </span><span texencoding="V^k_j" style="vertical-align:-8px"><img src="" width="18.5" height="22" /></span><span>, </span><span texencoding="LB^k_j" style="vertical-align:-8px"><img src="" width="27.5" height="22" /></span><span> and </span><span texencoding="UB^k_j" style="vertical-align:-8px"><img src="" width="28.5" height="22" /></span><span> are respectively the flux (in mmol/h), lower bound (in mmol/h/gdw) and upper bound (in mmol/h/gdw) for reaction </span><span style=' font-style: italic;'>j</span><span> for organism </span><span style=' font-style: italic;'>k</span><span>, </span><span texencoding="u^{com}_i" style="vertical-align:-8px"><img src="" width="27" height="22" /></span><span> is the community uptake bound for metabolite </span><span style=' font-style: italic;'>i</span><span>, </span><span texencoding="X^k" style="vertical-align:-5px"><img src="" width="18.5" height="19" /></span><span> is the biomass (in gdw) of organism </span><span style=' font-style: italic;'>k</span><span>, </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">μ</span><span> is the community growth rate, </span><span texencoding="\textbf{I}^k" style="vertical-align:-5px"><img src="" width="13.5" height="19" /></span><span> is the set of metabolites of organism </span><span style=' font-style: italic;'>k</span><span>, </span><span texencoding="\textbf{I}^{com}" style="vertical-align:-5px"><img src="" width="25.5" height="19" /></span><span> is the set of community metabolites in the community exchange space, </span><span texencoding="\textbf{J}^k" style="vertical-align:-5px"><img src="" width="15" height="19" /></span><span> is the set of reactions for organism k, </span><span texencoding="\textbf{K}" style="vertical-align:-5px"><img src="" width="14.5" height="18" /></span><span> is the set of organisms in the community, and </span><span texencoding="ex(i) \in \textbf{J}^k" style="vertical-align:-5px"><img src="" width="63" height="20" /></span><span> is the exchange reaction in organism </span><span style=' font-style: italic;'>k</span><span> for extracellular metabolite </span><span style=' font-style: italic;'>i</span><span>. See ref. [1] for the derivation and detailed explanation.</span></div><div  class = 'S2'><span>Throughout the tutorial, using a hypothetical model of four </span><span style=' font-style: italic;'>E. coli</span><span> mutants auxotrophic for amino acids, we will demonstrate the three different functionalities of the module: (1) computing the maximum community growth rate using the function SteadyCom.m, (2) performing flux variability analysis under a given community growth rate using SteadyComFVA.m, and (3) analyzing the pairwise relationship between flux/biomass variables using a technique similar to Pareto-optimal analysis by calling the function SteadyComPOA.m</span></div><h2  class = 'S1'><span>EQUIPMENT SETUP</span></h2><div  class = 'S2'><span>If necessary, initialise the cobra toolbox and select a solver by running:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S4'><span style="white-space: pre"><span >initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="D2D3B12E" data-testid="output_0" data-width="420" data-height="899" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">      _____   _____   _____   _____     _____     |
     /  ___| /  _  \ |  _  \ |  _  \   / ___ \    |   COnstraint-Based Reconstruction and Analysis
     | |     | | | | | |_| | | |_| |  | |___| |   |   The COBRA Toolbox - 2017
     | |     | | | | |  _  { |  _  /  |  ___  |   |
     | |___  | |_| | | |_| | | | \ \  | |   | |   |   Documentation:
     \_____| \_____/ |_____/ |_|  \_\ |_|   |_|   |   <a href="http://opencobra.github.io/cobratoolbox" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">http://opencobra.github.io/cobratoolbox</a>
                                                  | 

 &gt; Checking if git is installed ...  Done.
 &gt; Checking if the repository is tracked using git ...  Done.
 &gt; Checking if curl is installed ...  Done.
 &gt; Checking if remote can be reached ...  Done.
 &gt; Initializing and updating submodules ... Done.
 &gt; Adding all the files of The COBRA Toolbox ...  Done.
 &gt; Define CB map output... set to svg.
 &gt; Retrieving models ...   Done.
 &gt; TranslateSBML is installed and working properly.
 &gt; Configuring solver environment variables ...
   - [*---] ILOG_CPLEX_PATH: /Applications/IBM/ILOG/CPLEX_Studio1271/cplex/matlab/x86-64_osx
   - [*---] GUROBI_PATH: /Library/gurobi650/mac64/matlab
   - [----] TOMLAB_PATH :  --&gt; set this path manually after installing the solver ( see <a href="https://opencobra.github.io/cobratoolbox/docs/solvers.html" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">instructions</a> )
   - [----] MOSEK_PATH :  --&gt; set this path manually after installing the solver ( see <a href="https://opencobra.github.io/cobratoolbox/docs/solvers.html" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">instructions</a> )
   Done.
 &gt; Checking available solvers and solver interfaces ... Done.
 &gt; Setting default solvers ... Done.
 &gt; Saving the MATLAB path ... Done.
   - The MATLAB path was saved in the default location.

 &gt; Summary of available solvers and solver interfaces

			Support 	   LP 	 MILP 	   QP 	 MIQP 	  NLP
	----------------------------------------------------------------------
	cplex_direct 	active        	    0 	    0 	    0 	    0 	    -
	dqqMinos     	active        	    1 	    - 	    - 	    - 	    -
	glpk         	active        	    1 	    1 	    - 	    - 	    -
	gurobi       	active        	    1 	    1 	    1 	    1 	    -
	ibm_cplex    	active        	    1 	    1 	    1 	    - 	    -
	matlab       	active        	    1 	    - 	    - 	    - 	    1
	mosek        	active        	    0 	    0 	    0 	    - 	    -
	pdco         	active        	    1 	    - 	    1 	    - 	    -
	quadMinos    	active        	    1 	    - 	    - 	    - 	    1
	tomlab_cplex 	active        	    0 	    0 	    0 	    0 	    -
	qpng         	passive       	    - 	    - 	    1 	    - 	    -
	tomlab_snopt 	passive       	    - 	    - 	    - 	    - 	    0
	gurobi_mex   	legacy        	    0 	    0 	    0 	    0 	    -
	lindo_old    	legacy        	    0 	    - 	    - 	    - 	    -
	lindo_legacy 	legacy        	    0 	    - 	    - 	    - 	    -
	lp_solve     	legacy        	    1 	    - 	    - 	    - 	    -
	opti         	legacy        	    0 	    0 	    0 	    0 	    0
	----------------------------------------------------------------------
	Total        	-             	    8 	    3 	    4 	    1 	    2

 + Legend: - = not applicable, 0 = solver not compatible or not installed, 1 = solver installed.


 &gt; You can solve LP problems using: 'dqqMinos' - 'glpk' - 'gurobi' - 'ibm_cplex' - 'matlab' - 'pdco' - 'quadMinos' - 'lp_solve' 
 &gt; You can solve MILP problems using: 'glpk' - 'gurobi' - 'ibm_cplex' 
 &gt; You can solve QP problems using: 'gurobi' - 'ibm_cplex' - 'pdco' - 'qpng' 
 &gt; You can solve MIQP problems using: 'gurobi' 
 &gt; You can solve NLP problems using: 'matlab' - 'quadMinos' 

 &gt; Checking for available updates ...
 --&gt; You cannot update your fork using updateCobraToolbox(). [97ac46 @ master].
     Please use the MATLAB.devTools (<a href="https://github.com/opencobra/MATLAB.devTools" style="white-space: pre; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">https://github.com/opencobra/MATLAB.devTools</a>).</div></div></div></div></div><div  class = 'S6'><span>All SteadyCom functions involve only solving linear programming problems. Any solvers supported by the COBRA toolbox will work. But SteadyCom contains specialized codes for IBM ILOG Cplex which was tested to run significantly faster for SteadyComFVA and SteadyComPOA for larger problems through calling the Cplex object in Matlab directly. </span></div><div  class = 'S2'><span>Please note that parallelization requires a working installation of the Parallel Computing Toolbox.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S4'><span style="white-space: pre"><span >changeCobraSolver(</span><span style="color: rgb(170, 4, 249);">'ibm_cplex'</span><span >, </span><span style="color: rgb(170, 4, 249);">'LP'</span><span >);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="39CE82B4" data-testid="output_1" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> &gt; IBM ILOG CPLEX interface added to MATLAB path.</div></div></div></div></div><h2  class = 'S7'><span>PROCEDURE</span></h2><h2  class = 'S1'><span>Model Construction</span></h2><div  class = 'S2'><span>Load the </span><span style=' font-style: italic;'>E. coli</span><span> iAF1260 model in the COBRA toolbox.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">global </span><span >CBTDIR</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >iAF1260 = readCbModel([CBTDIR filesep </span><span style="color: rgb(170, 4, 249);">'test' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'models' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'iAF1260.mat'</span><span >]);</span></span></div></div></div><div  class = 'S6'><span>Polish the model a little bit:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% convert the compartment format from e.g., '_c' to '[c]'</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >iAF1260.mets = regexprep(iAF1260.mets, </span><span style="color: rgb(170, 4, 249);">'_([^_]+)$'</span><span >, </span><span style="color: rgb(170, 4, 249);">'\[$1\]'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% make all empty cells in cell arrays to be empty string</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >fieldToBeCellStr = {</span><span style="color: rgb(170, 4, 249);">'metFormulas'</span><span >; </span><span style="color: rgb(170, 4, 249);">'genes'</span><span >; </span><span style="color: rgb(170, 4, 249);">'grRules'</span><span >; </span><span style="color: rgb(170, 4, 249);">'metNames'</span><span >; </span><span style="color: rgb(170, 4, 249);">'rxnNames'</span><span >; </span><span style="color: rgb(170, 4, 249);">'subSystems'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >j = 1:numel(fieldToBeCellStr)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    iAF1260.(fieldToBeCellStr{j})(cellfun(@isempty, iAF1260.(fieldToBeCellStr{j}))) = {</span><span style="color: rgb(170, 4, 249);">''</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S6'><span>Add a methionine export reaction to allow the export of methionine.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S4'><span style="white-space: pre"><span >iAF1260 = addReaction(iAF1260,{</span><span style="color: rgb(170, 4, 249);">'METt3pp'</span><span >,</span><span style="color: rgb(170, 4, 249);">''</span><span >},</span><span style="color: rgb(170, 4, 249);">'met__L[c] + h[c] =&gt; met__L[p] + h[p]'</span><span >);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A5D57625" data-testid="output_2" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">METt3pp	h[c] + met__L[c] 	-&gt;	h[p] + met__L[p] </div></div></div></div></div><div  class = 'S6'><span>Reactions essential for amino acid autotrophy:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >argH = {</span><span style="color: rgb(170, 4, 249);">'ARGSL'</span><span >};  </span><span style="color: rgb(2, 128, 9);">% essential for arginine biosynthesis</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >lysA = {</span><span style="color: rgb(170, 4, 249);">'DAPDC'</span><span >};  </span><span style="color: rgb(2, 128, 9);">% essential for lysine biosynthesis</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >metA = {</span><span style="color: rgb(170, 4, 249);">'HSST'</span><span >};  </span><span style="color: rgb(2, 128, 9);">% essential for methionine biosynthesis</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >ilvE = {</span><span style="color: rgb(170, 4, 249);">'PPNDH'</span><span >};  </span><span style="color: rgb(2, 128, 9);">% essential for phenylalanine biosynthesis</span></span></div></div></div><div  class = 'S6'><span>Reactions essential for exporting amino acids:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >argO = {</span><span style="color: rgb(170, 4, 249);">'ARGt3pp'</span><span >};  </span><span style="color: rgb(2, 128, 9);">% Evidence for an arginine exporter encoded by yggA (argO) that is regulated by the LysR-type transcriptional regulator ArgP in Escherichia coli.</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >lysO = {</span><span style="color: rgb(170, 4, 249);">'LYSt3pp'</span><span >};  </span><span style="color: rgb(2, 128, 9);">% Distinct paths for basic amino acid export in Escherichia coli: YbjE (LysO) mediates export of L-lysine</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >yjeH = {</span><span style="color: rgb(170, 4, 249);">'METt3pp'</span><span >};  </span><span style="color: rgb(2, 128, 9);">% YjeH is a novel L-methionine and branched chain amino acids exporter in Escherichia coli</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >yddG = {</span><span style="color: rgb(170, 4, 249);">'PHEt2rpp'</span><span >};  </span><span style="color: rgb(2, 128, 9);">% YddG from Escherichia coli promotes export of aromatic amino acids.</span></span></div></div></div><div  class = 'S6'><span>Now make four copies of the model with auxotrophy for different amino acids and inability to export amino acids:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% auxotrophic for Lys and Met, not exporting Phe</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >Ec1 = iAF1260;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >Ec1 = changeRxnBounds(Ec1, [lysA; metA; yddG], 0, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% auxotrophic for Arg and Phe, not exporting Met</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >Ec2 = iAF1260;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >Ec2 = changeRxnBounds(Ec2, [argH; yjeH; ilvE], 0, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Auxotrophic for Arg and Phe, not exporting Lys</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >Ec3 = iAF1260;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >Ec3 = changeRxnBounds(Ec3, [argH; lysO; ilvE], 0, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Auxotrophic for Lys and Met, not exporting Arg</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >Ec4 = iAF1260;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >Ec4 = changeRxnBounds(Ec4, [argO; lysA; metA], 0, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div></div><div  class = 'S6'><span>Now none of the four organisms can grow alone and they must cross feed each other to survive. See Figure 1 in ref. </span><a href = "http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005539"><span>[1]</span></a><span> for the visualization of the community. </span></div><div  class = 'S2'><span>Get the extracellular metabolites, the corresponding exchange reactions and the uptake rates for the </span><span style=' font-style: italic;'>E. coli</span><span> model, which are used later to constrain the community model:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% extracellular metabolites (met[e])</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >metEx = strcmp(getCompartment(iAF1260.mets),</span><span style="color: rgb(170, 4, 249);">'e'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% the corresponding exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >rxnExAll = find(sum(iAF1260.S ~= 0, 1) == 1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >[rxnEx, ~] = find(iAF1260.S(metEx, rxnExAll)');  </span><span style="color: rgb(2, 128, 9);">% need to be in the same order as metEx</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >rxnEx = rxnExAll(rxnEx);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% exchange rate</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >lbEx = iAF1260.lb(rxnEx);</span></span></div></div></div><div  class = 'S6'><span>Create a community model with the four </span><span style=' font-style: italic;'>E. coli</span><span> tagged as 'Ec1', 'Ec2', 'Ec3', 'Ec4' respectively by calling </span><span style=' font-family: monospace;'>createMultipleSpeciesModel</span><span>. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >nameTagsModel = {</span><span style="color: rgb(170, 4, 249);">'Ec1'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec2'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec3'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec4'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = createMultipleSpeciesModel({Ec1; Ec2; Ec3; Ec4}, nameTagsModel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom.csense = char(</span><span style="color: rgb(170, 4, 249);">'E' </span><span >* ones(1,numel(EcCom.mets)));  </span><span style="color: rgb(2, 128, 9);">% correct the csense</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >clear </span><span style="color: rgb(170, 4, 249);">Ec1 Ec2 Ec3 Ec4</span></span></div></div></div><div  class = 'S6'><span>The model </span><span style=' font-family: monospace;'>EcCom</span><span> contains a community compartment denoted by </span><span style=' font-family: monospace;'>[u]</span><span> to allow exchange between organisms. Each organism-specific reaction/metabolite is prepended with the corresponding tag.</span></div><div  class = 'S2'><span>Retreive the names and ids for organism/community exchange reactions/metabolites which are necessary for computation:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[EcCom.infoCom, EcCom.indCom] = getMultiSpeciesModelId(EcCom, nameTagsModel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >disp(EcCom.infoCom);</span></span></div></div></div><div  class = 'S6'><span style=' font-family: monospace;'>EcCom.infoCom </span><span>contains reaction/metabolite names (from </span><span style=' font-family: monospace;'>EcCom.rxns</span><span>/</span><span style=' font-family: monospace;'>EcCom.mets</span><span>) for the community exchange reactions (</span><span style=' font-family: monospace;'>*.EXcom</span><span>), organism-community exchange reactions (</span><span style=' font-family: monospace;'>*.EXsp</span><span>), community metabolites (</span><span style=' font-family: monospace;'>*.Mcom</span><span>), organism-specific extracellular metabolite (</span><span style=' font-family: monospace;'>*.Msp</span><span>). If a host model is specified, there will also be non-empty </span><span style=' font-family: monospace;'>*.EXhost</span><span> and </span><span style=' font-family: monospace;'>*.Mhost </span><span>for the host-specific exchange reactions and metabolites. The fields </span><span style=' font-family: monospace;'>*.rxnSps</span><span>/</span><span style=' font-family: monospace;'>*.metSps</span><span> give information on which organism a reaction/metabolite belongs to.</span></div><div  class = 'S2'><span style=' font-family: monospace;'>indCom </span><span>has the same structure as </span><span style=' font-family: monospace;'>infoCom</span><span> but contains the indices rather than names. </span><span style=' font-family: monospace;'>infoCom</span><span> and </span><span style=' font-family: monospace;'>indCom</span><span> are attached as fields of the model </span><span style=' font-family: monospace;'>EcCom</span><span> because SteadyCom requires this information from the input model for computation. Incorporate also the names and indices for the biomass reactions which are necessary for computing growth:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >rxnBiomass = strcat(nameTagsModel, </span><span style="color: rgb(170, 4, 249);">'BIOMASS_Ec_iAF1260_core_59p81M'</span><span >);  </span><span style="color: rgb(2, 128, 9);">% biomass reaction names</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >rxnBiomassId = findRxnIDs(EcCom, rxnBiomass);  </span><span style="color: rgb(2, 128, 9);">% ids</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom.infoCom.spBm = rxnBiomass;  </span><span style="color: rgb(2, 128, 9);">% .spBm for organism biomass reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >EcCom.indCom.spBm = rxnBiomassId;</span></span></div></div></div><div  class = 'S6'><span></span></div><h2  class = 'S7'><span>Finding Maximum Growth Rate Using SteadyCom</span></h2><div  class = 'S2'><span>Set community and organism-specific uptake rates to be the same as in the orginal iAF1260 model:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >[yn, id] = ismember(strrep(iAF1260.mets(metEx), </span><span style="color: rgb(170, 4, 249);">'[e]'</span><span >, </span><span style="color: rgb(170, 4, 249);">'[u]'</span><span >), EcCom.infoCom.Mcom);  </span><span style="color: rgb(2, 128, 9);">% map the metabolite name</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >assert(all(yn));  </span><span style="color: rgb(2, 128, 9);">% must be a 1-to-1 mapping</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom.lb(EcCom.indCom.EXcom(:,1)) = lbEx(id);  </span><span style="color: rgb(2, 128, 9);">% assign community uptake bounds</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom.ub(EcCom.indCom.EXcom(:,1)) = 1e5;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >EcCom.lb(EcCom.indCom.EXsp) = repmat(lbEx(id), 1, 4);  </span><span style="color: rgb(2, 128, 9);">% assign organism-specific uptake bounds</span></span></div></div></div><div  class = 'S6'><span>Set maximum allowed organism-specific uptake rates for the cross-feeding amino acids:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% only allow to take up the amino acids that one is auxotrophic for</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >exRate = 1;  </span><span style="color: rgb(2, 128, 9);">% maximum uptake rate for cross feeding AAs</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Ec1</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = changeRxnBounds(EcCom, {</span><span style="color: rgb(170, 4, 249);">'Ec1IEX_arg__L[u]tr'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec1IEX_phe__L[u]tr'</span><span >}, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = changeRxnBounds(EcCom, {</span><span style="color: rgb(170, 4, 249);">'Ec1IEX_met__L[u]tr'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec1IEX_lys__L[u]tr'</span><span >}, -exRate, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Ec2</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = changeRxnBounds(EcCom, {</span><span style="color: rgb(170, 4, 249);">'Ec2IEX_arg__L[u]tr'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec2IEX_phe__L[u]tr'</span><span >}, -exRate, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = changeRxnBounds(EcCom, {</span><span style="color: rgb(170, 4, 249);">'Ec2IEX_met__L[u]tr'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec2IEX_lys__L[u]tr'</span><span >}, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Ec3</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = changeRxnBounds(EcCom, {</span><span style="color: rgb(170, 4, 249);">'Ec3IEX_arg__L[u]tr'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec3IEX_phe__L[u]tr'</span><span >}, -exRate, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = changeRxnBounds(EcCom, {</span><span style="color: rgb(170, 4, 249);">'Ec3IEX_met__L[u]tr'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec3IEX_lys__L[u]tr'</span><span >}, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Ec4</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = changeRxnBounds(EcCom, {</span><span style="color: rgb(170, 4, 249);">'Ec4IEX_arg__L[u]tr'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec4IEX_phe__L[u]tr'</span><span >}, 0, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom = changeRxnBounds(EcCom, {</span><span style="color: rgb(170, 4, 249);">'Ec4IEX_met__L[u]tr'</span><span >; </span><span style="color: rgb(170, 4, 249);">'Ec4IEX_lys__L[u]tr'</span><span >}, -exRate, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% allow production of anything for each member</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre"><span >EcCom.ub(EcCom.indCom.EXsp(:)) = 1000;</span></span></div></div></div><div  class = 'S6'><span>Before the calculation, print the community uptake bounds for checking using </span><span style=' font-family: monospace;'>printUptakeBoundCom</span><span>:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S4'><span style="white-space: pre"><span >printUptakeBoundCom(EcCom, 1);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="5DFAF1D7" data-testid="output_3" data-width="420" data-height="381" data-hashorizontaloverflow="true" style="width: 450px; max-height: 392px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">         Mets	Comm.     Ec1       Ec2       Ec3       Ec4       
 ( 53) arg__L	0         0         -1        -1        0         
    ( 60) ca2	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
   ( 62) cbl1	0.01      -0.01     -0.01     -0.01     -0.01     
     ( 67) cl	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    ( 69) co2	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
( 70) cobalt2	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    ( 76) cu2	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    (108) fe2	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    (109) fe3	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
 (144) glc__D	8         -8        -8        -8        -8        
    (167) h2o	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
      (169) h	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
      (186) k	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
 (194) lys__L	0         -1        0         0         -1        
 (208) met__L	0         -1        0         0         -1        
    (211) mg2	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    (214) mn2	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
   (216) mobd	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    (219) na1	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    (221) nh4	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
     (228) o2	18.5      -18.5     -18.5     -18.5     -18.5     
 (237) phe__L	0         0         -1        -1        0         
     (239) pi	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    (260) so4	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
  (280) tungs	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    
    (299) zn2	1e+06     -1e+06    -1e+06    -1e+06    -1e+06    </div></div></div></div></div><div  class = 'S6'><span>Values under 'Comm.' are the community uptake bounds (+ve for uptake) and values under 'Ec1' are the Ec1-specific uptake bounds (-ve for uptake). </span></div><div  class = 'S2'><span>Create an option structure for calling SteadyCom and call the function. There are a range of options available, including setting algorithmic parameters, fixing growth rates for members, adding additional linear constraints in a general format, e.g., for molecular crowding effect. See </span><span style=' font-family: monospace;'>help SteadyCom </span><span>for more options.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >options = struct();</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.GRguess = 0.5;  </span><span style="color: rgb(2, 128, 9);">% initial guess for max. growth rate</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.GRtol = 1e-6;  </span><span style="color: rgb(2, 128, 9);">% tolerance for final growth rate</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.algorithm = 1;  </span><span style="color: rgb(2, 128, 9);">% use the default algorithm (simple guessing for bounds, followed by matlab fzero)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >[sol, result] = SteadyCom(EcCom, options);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="BABAB0ED" data-testid="output_4" data-width="420" data-height="227" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Find maximum community growth rate..
Model feasible at maintenance. Time elapsed: 1 / 1 sec
Iter        LB   To test        UB  Time elapsed (iteration/total)
   1  0.000000  0.500000       Inf  0 / 1 sec
   2  0.500000  0.721279       Inf  4 / 5 sec
   3  0.721279  0.735372       Inf  0 / 5 sec
   4  0.735372  0.742726       Inf  0 / 5 sec
 
 Func-count    x          f(x)             Procedure
    2        0.735372  -0.000807615        initial
    3        0.735378   -0.00079987        interpolation
    4         0.73599  -1.26127e-06        interpolation
    5         0.73599  -1.26127e-06        interpolation
 
Zero found in the interval [0.735372, 0.742726]
Maximum community growth rate: 0.735990 (abs. error &lt; 1e-06).	Time elapsed: 21 sec</div></div></div></div></div><div  class = 'S6'><span>The algorithm is an iterative procedure to find the maximum biomass at a given growth rate and to determine the maximum growth rate that is feasible for the required total biomass (default 1 gdw). Here the algorithm used is the simple guessing for find upper and lower bounds (Iter 1 to 4 in the output) followed by Matlab </span><span style=' font-family: monospace;'>fzero</span><span> (starting from the line '</span><span style=' font-family: monospace;'>Func-count</span><span>') to locate the root. The maximum growth rate calculated is 0.73599 /h, stored in </span><span style=' font-family: monospace;'>result.GRmax</span><span>.</span><span style=' font-style: italic;'> </span></div><div  class = 'S2'><span>The biomass for each organism (in gdw) is given by</span><span style=' font-style: italic;'> </span><span style=' font-family: monospace;'>result.BM</span><span>:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >jSp = 1:4</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'X_%s:  %.6f\n'</span><span >, EcCom.infoCom.spAbbr{jSp}, result.BM(jSp));</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="773BE0CE" data-testid="output_5" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">X_Ec1:  0.253294</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A11F47A7" data-testid="output_6" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">X_Ec2:  0.324611</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="0531E544" data-testid="output_7" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">X_Ec3:  0.185004</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D82C44C8" data-testid="output_8" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">X_Ec4:  0.237093</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span >disp(result);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="52996B01" data-testid="output_9" data-width="420" data-height="129" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    GRmax: 0.7360
      vBM: [4×1 double]
       BM: [4×1 double]
       Ut: [299×1 double]
       Ex: [299×1 double]
     flux: [9831×1 double]
    iter0: [0 11.4198 0 9.9476e-14]
     iter: [4×6 double]
     stat: 'optimal'</div></div></div></div></div><div  class = 'S6'><span style=' font-family: monospace;'>result.vBM</span><span> contains the biomass production rates (in gdw / h), equal to </span><span style=' font-family: monospace;'>result.BM * result.GRmax. </span><span>Since the total community biomass is defaulted to be 1 gdw, the biomass for each organism coincides with its relative abundance. Note that the community uptake bounds in this sense are normalized per gdw of the community biomass. So the lower bound for the exchange reaction </span><span style=' font-family: monospace;'>EX_glc__D[u]</span><span> being 8 can be interpreted as the maximum amount of glucose available to the community being at a rate of 8 mmol per hour for 1 gdw of community biomass. Similarly, all fluxes in </span><span style=' font-family: monospace;'>result.flux </span><span>(</span><span texencoding="V^k_j" style="vertical-align:-8px"><img src="" width="18.5" height="22" /></span><span>)</span><span style=' font-family: monospace;'> </span><span>has the unit mmol / h / [gdw of comm. biomass]. It differs from the specific rate (traditionally denoted by </span><span texencoding="v^k_j" style="vertical-align:-8px"><img src="" width="14.5" height="22" /></span><span>) of an organism in the usual sense (in the unit of mmol / h / [gdw of organism biomass]) by </span><span texencoding="V^k_j=X^kv^k_j" style="vertical-align:-8px"><img src="" width="64" height="23" /></span><span> where </span><span texencoding="X^k" style="vertical-align:-5px"><img src="" width="18.5" height="19" /></span><span> is the biomass of the organism. </span><span style=' font-family: monospace;'>result.Ut</span><span style=' font-style: italic;'> </span><span>and</span><span style=' font-style: italic;'> </span><span style=' font-family: monospace;'>result.Ex </span><span>are the community uptake and export rates respectively, corresponding to the exchange reactions in </span><span style=' font-family: monospace;'>EcCom.infoCom.EXcom</span><span>.</span><span style=' font-family: monospace;'> </span></div><div  class = 'S2'><span style=' font-family: monospace;'>result.iter0 </span><span>is the info for solving the model at zero growth rate and </span><span style=' font-family: monospace;'>result.iter </span><span>records the info during iteration of the algorithm:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >iter = [0, result.iter0, NaN; result.iter];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >j = 0 : size(iter, 1)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >j == 0</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'#iter\tgrowth rate (mu)\tmax. biomass (sum(X))\tmu * sum(X)\tmax. infeasibility\tguess method\n'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >        fprintf(</span><span style="color: rgb(170, 4, 249);">'%5d\t%16.6f\t%21.6f\t%11.6f\t%18.6e\t%d\n'</span><span >, iter(j,:))</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="FBB81560" data-testid="output_10" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">#iter	growth rate (mu)	max. biomass (sum(X))	mu * sum(X)	max. infeasibility	guess method</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="B7746AB8" data-testid="output_11" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    0	        0.000000	            11.419845	   0.000000	      9.947598e-14	NaN</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="0DA8B635" data-testid="output_12" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    1	        0.500000	             1.442559	   0.721279	      3.493989e-10	0</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="D9D51985" data-testid="output_13" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    2	        0.721279	             1.019539	   0.735372	      3.668634e-10	0</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="78A2D588" data-testid="output_14" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    3	        0.735372	             1.000808	   0.735966	      1.706138e-10	0</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="4AB14698" data-testid="output_15" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    4	        0.742726	             0.000000	   0.000000	      0.000000e+00	2</div></div></div></div></div><div  class = 'S6'><span style=' font-family: monospace;'>mu * sum(X)</span><span> in the forth column is equal to the biomass production rate. </span></div><div  class = 'S2'><span>The fifth column contains the maximum infeasibility of the solutions in each iteration.</span></div><div  class = 'S2'><span>Guess method in the last column represents the method used for guessing the growth rate solved in the current iteration:</span></div><div  class = 'S2'><span>0: the default simple guess by </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi&gt;&amp;mu;&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;next&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi&gt;&amp;mu;&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;current&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;mrow&gt;&lt;msubsup&gt;&lt;mrow&gt;&lt;mo&gt;&amp;sum;&lt;/mo&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;k&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;K&lt;/mi&gt;&lt;/mrow&gt;&lt;/msubsup&gt;&lt;mrow&gt;&lt;msubsup&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;X&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;k&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;current&lt;/mi&gt;&lt;/mrow&gt;&lt;/msubsup&gt;&lt;/mrow&gt;&lt;/mrow&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-8px"><img src="" width="160" height="27" /></span><span> (</span><span style=' font-style: italic;'>K</span><span> is the total number of organisms)</span></div><div  class = 'S2'><span>1: bisection method</span></div><div  class = 'S2'><span>2: bisection or at least 1% away from the bounds if the simple guess is too close to the bounds (&lt;1%)</span></div><div  class = 'S2'><span>3. 1% away from the current growth rate if the simple guess is too close to the current growth rate</span></div><div  class = 'S2'><span>From the table, we can see that at the growth rate 0.742726 (iter 4), the max. biomass is 0, while at growth rate 0.735372, max. biomass = 1.0008 &gt; 1. Therefore we have both an lower and upper bound for the max. growth rate. Then fzero is initiated to solve for the max. growth rate that gives max. biomass &gt;= 1.</span></div><div  class = 'S2'><span>Two other algorithms for the iterative procedure are also implemented: simple guessing only and the bisection method. Compare their results with simple guessing + matlab fzero run above:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >options.algorithm = 2;  </span><span style="color: rgb(2, 128, 9);">% use the simple guessing algorithm</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >[sol2, result2] = SteadyCom(EcCom, options);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="7422FEDC" data-testid="output_16" data-width="420" data-height="297" data-hashorizontaloverflow="true" style="width: 450px; max-height: 308px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Find maximum community growth rate..
Model feasible at maintenance. Time elapsed: 1 / 1 sec
Iter        LB   To test        UB  Time elapsed (iteration/total)
   1  0.000000  0.500000       Inf  0 / 1 sec
   2  0.500000  0.721279       Inf  4 / 5 sec
   3  0.721279  0.735372       Inf  0 / 5 sec
   4  0.735372  0.742726       Inf  0 / 5 sec
   5  0.735372  0.739049  0.742726  0 / 5 sec
   6  0.735372  0.737211  0.739049  0 / 5 sec
   7  0.735372  0.736291  0.737211  0 / 5 sec
   8  0.735372  0.735832  0.736291  0 / 6 sec
   9  0.735832  0.736062  0.736291  1 / 7 sec
  10  0.735832  0.735947  0.736062  0 / 7 sec
  11  0.735947  0.736004  0.736062  1 / 8 sec
  12  0.735947  0.735975  0.736004  0 / 8 sec
  13  0.735975  0.735990  0.736004  2 / 10 sec
  14  0.735990  0.735997  0.736004  0 / 10 sec
  15  0.735990  0.735993  0.735997  0 / 10 sec
  16  0.735990  0.735991  0.735993  0 / 11 sec
  17  0.735990  0.735991  0.735991  0 / 11 sec
Maximum community growth rate: 0.735991 (abs. error &lt; 1e-06).	Time elapsed: 14 sec</div></div></div></div><div class="inlineWrapper"><div  class = 'S13'><span style="white-space: pre"><span >options.algorithm = 3;  </span><span style="color: rgb(2, 128, 9);">% use the bisection algorithm</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >[sol3, result3] = SteadyCom(EcCom, options);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="3C8BF9D1" data-testid="output_17" data-width="420" data-height="353" data-hashorizontaloverflow="true" style="width: 450px; max-height: 364px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Find maximum community growth rate..
Model feasible at maintenance. Time elapsed: 0 / 0 sec
Iter        LB   To test        UB  Time elapsed (iteration/total)
   1  0.000000  0.500000       Inf  0 / 0 sec
   2  0.500000  1.000000       Inf  3 / 4 sec
   3  0.500000  0.750000  1.000000  0 / 4 sec
   4  0.500000  0.625000  0.750000  4 / 8 sec
   5  0.625000  0.687500  0.750000  5 / 13 sec
   6  0.687500  0.718750  0.750000  0 / 13 sec
   7  0.718750  0.734375  0.750000  0 / 13 sec
   8  0.734375  0.742188  0.750000  0 / 13 sec
   9  0.734375  0.738281  0.742188  0 / 13 sec
  10  0.734375  0.736328  0.738281  0 / 13 sec
  11  0.734375  0.735352  0.736328  0 / 14 sec
  12  0.735352  0.735840  0.736328  0 / 14 sec
  13  0.735840  0.736084  0.736328  0 / 14 sec
  14  0.735840  0.735962  0.736084  0 / 15 sec
  15  0.735962  0.736023  0.736084  1 / 16 sec
  16  0.735962  0.735992  0.736023  0 / 16 sec
  17  0.735962  0.735977  0.735992  0 / 17 sec
  18  0.735977  0.735985  0.735992  2 / 18 sec
  19  0.735985  0.735989  0.735992  0 / 19 sec
  20  0.735989  0.735991  0.735992  0 / 19 sec
  21  0.735991  0.735991  0.735992  0 / 19 sec
Maximum community growth rate: 0.735991 (abs. error &lt; 1e-06).	Time elapsed: 26 sec</div></div></div></div></div><div  class = 'S6'><span>The time used for each algorithm in the tested machine is:</span></div><div  class = 'S2'><span>(1) simple guess for bounds followed by Matlab fzero: 18 sec</span></div><div  class = 'S2'><span>(2) simple guess alone: 35 sec</span></div><div  class = 'S2'><span>(3) bisection: 70 sec</span></div><div  class = 'S2'><span>Algorithm (1) appears to be the fastest in most case although the simple guess algorithm can sometimes also outperform it. The most conservative bisection method can already guarantee convergence within around 20 iterations, i.e., solving ~20 LPs for an optimality gap (</span><span style=' font-family: monospace;'>options.GRtol</span><span>) of 1e-6.</span></div><div  class = 'S2'><span></span></div><h2  class = 'S7'><span>Analyzing Flux Variability Using SteadyComFVA</span></h2><div  class = 'S2'><span>Now we want to analyze the variability of the organism abundance at various growth rates. Choose more options and call </span><span style=' font-family: monospace;'>SteadyComFVA</span><span>:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% percentage of maximum total biomass of the community required. 100 for sum(biomass) = 1 (1 is the default total biomass)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.optBMpercent = 100;  </span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >n = size(EcCom.S, 2);  </span><span style="color: rgb(2, 128, 9);">% number of reactions in the model</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% options.rxnNameList is the list of reactions subject to FVA. Can be reaction names or indices.</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Use n + j for the biomass variable of the j-th organism. Alternatively, use {'X_j'} </span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% for biomass variable of the j-th organism or {'X_Ec1'} for Ec1 (the abbreviation in EcCom.infoCom.spAbbr)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.rxnNameList = {</span><span style="color: rgb(170, 4, 249);">'X_Ec1'</span><span >; </span><span style="color: rgb(170, 4, 249);">'X_Ec2'</span><span >; </span><span style="color: rgb(170, 4, 249);">'X_Ec3'</span><span >; </span><span style="color: rgb(170, 4, 249);">'X_Ec4'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.optGRpercent = [89:0.2:99, 99.1:0.1:100];  </span><span style="color: rgb(2, 128, 9);">% perform FVA at various percentages of the maximum growth rate, 89, 89.1, 89.2, ..., 100</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >[fvaComMin,fvaComMax] = SteadyComFVA(EcCom, options);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="C80EB6D9" data-testid="output_18" data-width="420" data-height="591" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Find maximum community growth rate..
Model feasible at maintenance. Time elapsed: 1 / 1 sec
Iter        LB   To test        UB  Time elapsed (iteration/total)
   1  0.000000  0.500000       Inf  0 / 1 sec
   2  0.500000  1.000000       Inf  4 / 5 sec
   3  0.500000  0.750000  1.000000  0 / 5 sec
   4  0.500000  0.625000  0.750000  5 / 11 sec
   5  0.625000  0.687500  0.750000  7 / 17 sec
   6  0.687500  0.718750  0.750000  0 / 17 sec
   7  0.718750  0.734375  0.750000  0 / 17 sec
   8  0.734375  0.742188  0.750000  0 / 18 sec
   9  0.734375  0.738281  0.742188  0 / 18 sec
  10  0.734375  0.736328  0.738281  0 / 18 sec
  11  0.734375  0.735352  0.736328  0 / 18 sec
  12  0.735352  0.735840  0.736328  0 / 19 sec
  13  0.735840  0.736084  0.736328  0 / 19 sec
  14  0.735840  0.735962  0.736084  0 / 19 sec
  15  0.735962  0.736023  0.736084  2 / 21 sec
  16  0.735962  0.735992  0.736023  0 / 21 sec
  17  0.735962  0.735977  0.735992  0 / 22 sec
  18  0.735977  0.735985  0.735992  2 / 24 sec
  19  0.735985  0.735989  0.735992  0 / 24 sec
  20  0.735989  0.735991  0.735992  0 / 24 sec
  21  0.735991  0.735991  0.735992  0 / 24 sec
Maximum community growth rate: 0.735991 (abs. error &lt; 1e-06).	Time elapsed: 33 sec

FVA for 4 sets of fluxes/biomass at growth rate 0.655032 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.044053	 0.787578
   2	  50	     X_Ec2	 0.038253	 0.720492
   3	  75	     X_Ec3	 0.021200	 0.696956
   4	 100	     X_Ec4	 0.029222	 0.697238
BMmax adjustment: 1
BMmax adjustment: 2
BMmax adjustment: 3
BMmax adjustment: 4
BMmax adjustment: 5
BMmax adjustment: 6
BMmax adjustment: 7
BMmax adjustment: 8
BMmax adjustment: 9
BMmax adjustment: 10</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="93A02C33" data-testid="output_19" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: pre-wrap; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;">Warning: Model not feasible.</div><div class="diagnosticMessage-stackPart" style="white-space: pre; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="F5E64D27" data-testid="output_20" data-width="420" data-height="2621" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">FVA for 4 sets of fluxes/biomass at growth rate 0.657976 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.046186	 0.783368
   2	  50	     X_Ec2	 0.039919	 0.713899
   3	  75	     X_Ec3	 0.022092	 0.689206
   4	 100	     X_Ec4	 0.030498	 0.689833

FVA for 4 sets of fluxes/biomass at growth rate 0.659448 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.047304	 0.781210
   2	  50	     X_Ec2	 0.040788	 0.710505
   3	  75	     X_Ec3	 0.022556	 0.685205
   4	 100	     X_Ec4	 0.031163	 0.686016

FVA for 4 sets of fluxes/biomass at growth rate 0.660919 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.048458	 0.779016
   2	  50	     X_Ec2	 0.041682	 0.707043
   3	  75	     X_Ec3	 0.023033	 0.681117
   4	 100	     X_Ec4	 0.031848	 0.682120

FVA for 4 sets of fluxes/biomass at growth rate 0.662391 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.049649	 0.776783
   2	  50	     X_Ec2	 0.042603	 0.703511
   3	  75	     X_Ec3	 0.023523	 0.676937
   4	 100	     X_Ec4	 0.032553	 0.678142
BMmax adjustment: 1

FVA for 4 sets of fluxes/biomass at growth rate 0.663863 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.050880	 0.774509
   2	  50	     X_Ec2	 0.043552	 0.699897
   3	  75	     X_Ec3	 0.024028	 0.672653
   4	 100	     X_Ec4	 0.033283	 0.674078

FVA for 4 sets of fluxes/biomass at growth rate 0.665335 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.052152	 0.772192
   2	  50	     X_Ec2	 0.044530	 0.696203
   3	  75	     X_Ec3	 0.024547	 0.668265
   4	 100	     X_Ec4	 0.034036	 0.669928

FVA for 4 sets of fluxes/biomass at growth rate 0.666807 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.053466	 0.769834
   2	  50	     X_Ec2	 0.045538	 0.692431
   3	  75	     X_Ec3	 0.025082	 0.663776
   4	 100	     X_Ec4	 0.034812	 0.665686

FVA for 4 sets of fluxes/biomass at growth rate 0.668279 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.054825	      NaN
   2	  50	     X_Ec2	 0.046576	      NaN
   3	  75	     X_Ec3	      NaN	 0.659181
   4	 100	     X_Ec4	 0.035612	 0.661351

FVA for 4 sets of fluxes/biomass at growth rate 0.669751 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.056231	 0.764988
   2	  50	     X_Ec2	 0.047646	 0.684644
   3	  75	     X_Ec3	 0.026197	      NaN
   4	 100	     X_Ec4	 0.036437	 0.656920
BMmax adjustment: 1
BMmax adjustment: 2
BMmax adjustment: 3
BMmax adjustment: 4
BMmax adjustment: 5
BMmax adjustment: 6

FVA for 4 sets of fluxes/biomass at growth rate 0.671223 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.057686	      NaN
   2	  50	     X_Ec2	 0.048750	 0.680624
   3	  75	     X_Ec3	 0.026779	      NaN
   4	 100	     X_Ec4	 0.037288	 0.652387

FVA for 4 sets of fluxes/biomass at growth rate 0.672695 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.059191	 0.759959
   2	  50	     X_Ec2	 0.049888	 0.676516
   3	  75	     X_Ec3	 0.027379	      NaN
   4	 100	     X_Ec4	 0.038166	 0.647752

FVA for 4 sets of fluxes/biomass at growth rate 0.674167 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.060750	      NaN
   2	  50	     X_Ec2	 0.051063	 0.672316
   3	  75	     X_Ec3	 0.027996	      NaN
   4	 100	     X_Ec4	 0.039073	 0.643008

FVA for 4 sets of fluxes/biomass at growth rate 0.675639 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.062365	      NaN
   2	  50	     X_Ec2	 0.052275	 0.668022
   3	  75	     X_Ec3	 0.028632	 0.634496
   4	 100	     X_Ec4	 0.040009	      NaN

FVA for 4 sets of fluxes/biomass at growth rate 0.677111 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.064038	 0.752047
   2	  50	     X_Ec2	 0.053526	 0.663629
   3	  75	     X_Ec3	 0.029287	      NaN
   4	 100	     X_Ec4	 0.040976	 0.633183

FVA for 4 sets of fluxes/biomass at growth rate 0.678583 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.065772	 0.749305
   2	  50	     X_Ec2	 0.054818	 0.659135
   3	  75	     X_Ec3	 0.029963	 0.623739
   4	 100	     X_Ec4	 0.041975	 0.628092

FVA for 4 sets of fluxes/biomass at growth rate 0.680055 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.067571	 0.746507
   2	  50	     X_Ec2	 0.056153	 0.654536
   3	  75	     X_Ec3	 0.030659	 0.618150
   4	 100	     X_Ec4	 0.043007	 0.622877
BMmax adjustment: 1
BMmax adjustment: 2
BMmax adjustment: 3

FVA for 4 sets of fluxes/biomass at growth rate 0.681527 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.069437	      NaN
   2	  50	     X_Ec2	 0.057533	 0.649827
   3	  75	     X_Ec3	 0.031377	 0.612415
   4	 100	     X_Ec4	 0.044075	 0.617533

FVA for 4 sets of fluxes/biomass at growth rate 0.682999 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.071373	      NaN
   2	  50	     X_Ec2	 0.058959	 0.645006
   3	  75	     X_Ec3	 0.032118	 0.606527
   4	 100	     X_Ec4	 0.045179	 0.612055

FVA for 4 sets of fluxes/biomass at growth rate 0.684471 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.073384	      NaN
   2	  50	     X_Ec2	 0.060434	 0.640067
   3	  75	     X_Ec3	 0.032882	 0.600479
   4	 100	     X_Ec4	 0.046322	 0.606437

FVA for 4 sets of fluxes/biomass at growth rate 0.685943 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.075473	 0.734721
   2	  50	     X_Ec2	 0.061960	 0.635005
   3	  75	     X_Ec3	 0.033672	 0.594264
   4	 100	     X_Ec4	 0.047505	 0.600674

FVA for 4 sets of fluxes/biomass at growth rate 0.687415 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.077644	 0.731615
   2	  50	     X_Ec2	 0.063539	 0.629817
   3	  75	     X_Ec3	 0.034486	 0.587876
   4	 100	     X_Ec4	 0.048731	 0.594760

FVA for 4 sets of fluxes/biomass at growth rate 0.688887 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.079901	 0.728440
   2	  50	     X_Ec2	 0.065174	 0.624497
   3	  75	     X_Ec3	 0.035328	 0.581308
   4	 100	     X_Ec4	 0.050000	 0.588689

FVA for 4 sets of fluxes/biomass at growth rate 0.690359 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.082249	 0.725194
   2	  50	     X_Ec2	 0.066868	 0.619039
   3	  75	     X_Ec3	 0.036197	 0.574550
   4	 100	     X_Ec4	 0.051316	 0.582454

FVA for 4 sets of fluxes/biomass at growth rate 0.691831 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.084698	 0.721873
   2	  50	     X_Ec2	 0.068624	 0.613425
   3	  75	     X_Ec3	 0.037096	 0.567595
   4	 100	     X_Ec4	 0.052681	 0.576024
BMmax adjustment: 1
BMmax adjustment: 2
BMmax adjustment: 3
BMmax adjustment: 4
BMmax adjustment: 5
BMmax adjustment: 6
BMmax adjustment: 7
BMmax adjustment: 8
BMmax adjustment: 9
BMmax adjustment: 10</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="F98C3432" data-testid="output_21" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: pre-wrap; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;">Warning: Model not feasible.</div><div class="diagnosticMessage-stackPart" style="white-space: pre; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="B9F1F8BE" data-testid="output_22" data-width="420" data-height="143" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">BMmax adjustment: 1
BMmax adjustment: 2
BMmax adjustment: 3
BMmax adjustment: 4
BMmax adjustment: 5
BMmax adjustment: 6
BMmax adjustment: 7
BMmax adjustment: 8
BMmax adjustment: 9
BMmax adjustment: 10</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="857D9029" data-testid="output_23" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: pre-wrap; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;">Warning: Model not feasible.</div><div class="diagnosticMessage-stackPart" style="white-space: pre; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="A6F83DE1" data-testid="output_24" data-width="420" data-height="759" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">BMmax adjustment: 1
BMmax adjustment: 2
BMmax adjustment: 3
BMmax adjustment: 4
BMmax adjustment: 5
BMmax adjustment: 6
BMmax adjustment: 7
BMmax adjustment: 8
BMmax adjustment: 9

FVA for 4 sets of fluxes/biomass at growth rate 0.696247 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.092676	 0.711435
   2	  50	     X_Ec2	 0.074290	 0.595651
   3	  75	     X_Ec3	 0.039980	 0.545450
   4	 100	     X_Ec4	 0.057093	 0.555620

FVA for 4 sets of fluxes/biomass at growth rate 0.697719 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.095566	 0.707786
   2	  50	     X_Ec2	 0.076323	 0.589407
   3	  75	     X_Ec3	 0.041009	 0.537609
   4	 100	     X_Ec4	 0.058679	 0.548420

FVA for 4 sets of fluxes/biomass at growth rate 0.699191 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.098582	      NaN
   2	  50	     X_Ec2	 0.078435	 0.583010
   3	  75	     X_Ec3	 0.042075	 0.529518
   4	 100	     X_Ec4	 0.060328	 0.541006

FVA for 4 sets of fluxes/biomass at growth rate 0.700663 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.101732	 0.700210
   2	  50	     X_Ec2	 0.080630	 0.576441
   3	  75	     X_Ec3	 0.043179	 0.521166
   4	 100	     X_Ec4	 0.062043	 0.533368

FVA for 4 sets of fluxes/biomass at growth rate 0.702135 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.105023	 0.696276
   2	  50	     X_Ec2	 0.082912	 0.569710
   3	  75	     X_Ec3	 0.044323	 0.512540
   4	 100	     X_Ec4	 0.063828	 0.525494
BMmax adjustment: 1
BMmax adjustment: 2
BMmax adjustment: 3
BMmax adjustment: 4
BMmax adjustment: 5
BMmax adjustment: 6
BMmax adjustment: 7
BMmax adjustment: 8
BMmax adjustment: 9
BMmax adjustment: 10</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="7C71379E" data-testid="output_25" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: pre-wrap; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;">Warning: Model not feasible.</div><div class="diagnosticMessage-stackPart" style="white-space: pre; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="4DFD6154" data-testid="output_26" data-width="420" data-height="2803" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">BMmax adjustment: 1
BMmax adjustment: 2

FVA for 4 sets of fluxes/biomass at growth rate 0.705079 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.112067	      NaN
   2	  50	     X_Ec2	 0.087757	 0.555814
   3	  75	     X_Ec3	 0.046739	 0.494406
   4	 100	     X_Ec4	 0.067624	 0.508993

FVA for 4 sets of fluxes/biomass at growth rate 0.706551 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.115837	 0.683829
   2	  50	     X_Ec2	 0.090331	 0.548563
   3	  75	     X_Ec3	 0.048016	 0.484867
   4	 100	     X_Ec4	 0.069643	 0.500341

FVA for 4 sets of fluxes/biomass at growth rate 0.708023 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.119788	 0.679449
   2	  50	     X_Ec2	 0.093013	 0.541098
   3	  75	     X_Ec3	 0.049341	 0.474990
   4	 100	     X_Ec4	 0.071750	 0.491402

FVA for 4 sets of fluxes/biomass at growth rate 0.709495 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.123931	 0.674943
   2	  50	     X_Ec2	 0.095810	 0.533410
   3	  75	     X_Ec3	 0.050717	 0.464757
   4	 100	     X_Ec4	 0.073950	 0.482162

FVA for 4 sets of fluxes/biomass at growth rate 0.710967 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.128278	 0.670305
   2	  50	     X_Ec2	 0.098728	 0.525487
   3	  75	     X_Ec3	 0.052147	 0.454147
   4	 100	     X_Ec4	 0.076248	 0.472603

FVA for 4 sets of fluxes/biomass at growth rate 0.712439 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.132843	 0.665529
   2	  50	     X_Ec2	 0.101775	 0.517320
   3	  75	     X_Ec3	 0.053634	 0.443138
   4	 100	     X_Ec4	 0.078651	 0.462710
BMmax adjustment: 1
BMmax adjustment: 2
BMmax adjustment: 3
BMmax adjustment: 4
BMmax adjustment: 5
BMmax adjustment: 6

FVA for 4 sets of fluxes/biomass at growth rate 0.713911 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.137641	 0.660607
   2	  50	     X_Ec2	 0.104958	 0.508895
   3	  75	     X_Ec3	 0.055181	 0.431707
   4	 100	     X_Ec4	 0.081165	 0.452462

FVA for 4 sets of fluxes/biomass at growth rate 0.715383 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.142688	 0.655531
   2	  50	     X_Ec2	 0.108287	 0.500202
   3	  75	     X_Ec3	 0.056790	 0.419828
   4	 100	     X_Ec4	 0.083798	 0.441839

FVA for 4 sets of fluxes/biomass at growth rate 0.716855 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.148002	 0.650292
   2	  50	     X_Ec2	 0.111770	 0.491225
   3	  75	     X_Ec3	 0.058466	 0.407473
   4	 100	     X_Ec4	 0.086557	 0.430821

FVA for 4 sets of fluxes/biomass at growth rate 0.718327 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.153601	 0.644881
   2	  50	     X_Ec2	 0.115417	 0.481952
   3	  75	     X_Ec3	 0.060212	 0.394612
   4	 100	     X_Ec4	 0.089450	 0.419382
BMmax adjustment: 1

FVA for 4 sets of fluxes/biomass at growth rate 0.719799 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.159507	 0.639287
   2	  50	     X_Ec2	 0.119240	 0.472366
   3	  75	     X_Ec3	 0.062032	 0.381212
   4	 100	     X_Ec4	 0.092488	 0.407496

FVA for 4 sets of fluxes/biomass at growth rate 0.721271 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.165742	 0.633501
   2	  50	     X_Ec2	 0.123249	 0.462452
   3	  75	     X_Ec3	 0.063931	 0.367237
   4	 100	     X_Ec4	 0.095680	 0.395137

FVA for 4 sets of fluxes/biomass at growth rate 0.722743 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.172333	 0.627510
   2	  50	     X_Ec2	 0.127458	 0.452192
   3	  75	     X_Ec3	 0.065912	 0.352649
   4	 100	     X_Ec4	 0.099037	 0.382274
BMmax adjustment: 1
BMmax adjustment: 2

FVA for 4 sets of fluxes/biomass at growth rate 0.724215 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.179305	 0.621301
   2	  50	     X_Ec2	 0.131880	 0.441568
   3	  75	     X_Ec3	 0.067982	 0.337405
   4	 100	     X_Ec4	 0.102572	 0.368873

FVA for 4 sets of fluxes/biomass at growth rate 0.725687 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.186691	 0.614859
   2	  50	     X_Ec2	 0.136531	 0.430558
   3	  75	     X_Ec3	 0.070145	 0.321457
   4	 100	     X_Ec4	 0.106297	 0.354898

FVA for 4 sets of fluxes/biomass at growth rate 0.727159 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.194523	      NaN
   2	  50	     X_Ec2	 0.141428	 0.419142
   3	  75	     X_Ec3	 0.072407	 0.304754
   4	 100	     X_Ec4	 0.110228	 0.340309

FVA for 4 sets of fluxes/biomass at growth rate 0.728631 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.202839	 0.601215
   2	  50	     X_Ec2	 0.146588	 0.407296
   3	  75	     X_Ec3	 0.074774	 0.287239
   4	 100	     X_Ec4	 0.114380	 0.325063

FVA for 4 sets of fluxes/biomass at growth rate 0.729367 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.207190	 0.597632
   2	  50	     X_Ec2	 0.149273	 0.401204
   3	  75	     X_Ec3	 0.075999	 0.278158
   4	 100	     X_Ec4	 0.116544	 0.317179

FVA for 4 sets of fluxes/biomass at growth rate 0.730103 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.211679	 0.593976
   2	  50	     X_Ec2	 0.152032	 0.394995
   3	  75	     X_Ec3	 0.077253	 0.268849
   4	 100	     X_Ec4	 0.118771	 0.309112

FVA for 4 sets of fluxes/biomass at growth rate 0.730839 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.216310	 0.569878
   2	  50	     X_Ec2	 0.154868	 0.388666
   3	  75	     X_Ec3	 0.078538	 0.259305
   4	 100	     X_Ec4	 0.127080	 0.300856

FVA for 4 sets of fluxes/biomass at growth rate 0.731575 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.221090	 0.527616
   2	  50	     X_Ec2	 0.157783	 0.382212
   3	  75	     X_Ec3	 0.079852	 0.249515
   4	 100	     X_Ec4	 0.140974	 0.292403

FVA for 4 sets of fluxes/biomass at growth rate 0.732311 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.226026	 0.484427
   2	  50	     X_Ec2	 0.160780	 0.375631
   3	  75	     X_Ec3	 0.081199	 0.239469
   4	 100	     X_Ec4	 0.155428	 0.283745

FVA for 4 sets of fluxes/biomass at growth rate 0.733047 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.231124	 0.440276
   2	  50	     X_Ec2	 0.172784	 0.368917
   3	  75	     X_Ec3	 0.082578	 0.229158
   4	 100	     X_Ec4	 0.170469	 0.274876

FVA for 4 sets of fluxes/biomass at growth rate 0.733783 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.236391	 0.395127
   2	  50	     X_Ec2	 0.209556	 0.362068
   3	  75	     X_Ec3	 0.083992	 0.218570
   4	 100	     X_Ec4	 0.186124	 0.265787

FVA for 4 sets of fluxes/biomass at growth rate 0.734519 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.241835	 0.348944
   2	  50	     X_Ec2	 0.247095	 0.353601
   3	  75	     X_Ec3	 0.095040	 0.207693
   4	 100	     X_Ec4	 0.202424	 0.256468

FVA for 4 sets of fluxes/biomass at growth rate 0.735255 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.247466	 0.301686
   2	  50	     X_Ec2	 0.285430	 0.339473
   3	  75	     X_Ec3	 0.139450	 0.196515
   4	 100	     X_Ec4	 0.219401	 0.246911

FVA for 4 sets of fluxes/biomass at growth rate 0.735991 :
  No	   %	      Name	      Min	      Max
   1	  25	     X_Ec1	 0.253290	 0.253311
   2	  50	     X_Ec2	 0.324588	 0.324610
   3	  75	     X_Ec3	 0.185000	 0.185022
   4	 100	     X_Ec4	 0.237087	 0.237106</div></div></div></div></div><div  class = 'S6'><span>Similar to the output by </span><span style=' font-family: monospace;'>fluxVariability</span><span>, </span><span style=' font-family: monospace;'>fvaComMin</span><span> contains the minimum fluxes corresponding to the reactions in </span><span style=' font-family: monospace;'>options.rxnNameList</span><span>. </span><span style=' font-family: monospace;'>fvaComMax</span><span> contains the maximum fluxes. options.rxnNameList can be supplied as a (#rxns + #organism)-by-K matrix to analyze the variability of the K linear combinations of flux/biomass variables in the columns of the matrix. See </span><span style=' font-family: monospace;'>help SteadyComFVA</span><span> for more details.</span></div><div  class = 'S2'><span>We would also like to compare the results against the direct use of FBA and FVA by calling </span><span style=' font-family: monospace;'>optimizeCbModel</span><span> and </span><span style=' font-family: monospace;'>fluxVariability</span><span>:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >optGRpercentFBA = [89:2:99 99.1:0.1:100];  </span><span style="color: rgb(2, 128, 9);">% less dense interval to save time because the results are always the same for &lt; 99%</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >nGr = numel(optGRpercentFBA);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >[fvaFBAMin, fvaFBAMax] = deal(zeros(numel(options.rxnNameList), nGr));</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% change the objective function to the sum of all biomass reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom.c(:) = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom.c(EcCom.indCom.spBm) = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EcCom.csense = char(</span><span style="color: rgb(170, 4, 249);">'E' </span><span >* ones(1, numel(EcCom.mets)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >s = optimizeCbModel(EcCom);  </span><span style="color: rgb(2, 128, 9);">% run FBA</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >grFBA = s.f;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >jGr = 1:nGr</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'Growth rate %.4f :\n'</span><span >, grFBA * optGRpercentFBA(jGr)/100);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    [fvaFBAMin(:, jGr), fvaFBAMax(:, jGr)] = fluxVariability(EcCom, optGRpercentFBA(jGr), </span><span style="color: rgb(170, 4, 249);">'max'</span><span >, EcCom.infoCom.spBm, 2);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D67D6812" data-testid="output_27" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5091 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="25F22CDE" data-testid="output_28" data-width="420" data-height="31" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max
Starting parallel pool (parpool) using the 'local' profile ... connected to 2 workers.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D088EB60" data-testid="output_29" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5205 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D502613A" data-testid="output_30" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="35117DBA" data-testid="output_31" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5319 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="CC008DC8" data-testid="output_32" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="5DA8EE13" data-testid="output_33" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5434 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="BE6E2779" data-testid="output_34" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="735163D2" data-testid="output_35" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5548 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="E5B43A7C" data-testid="output_36" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="AF5AFD72" data-testid="output_37" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5663 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D9F2BF0A" data-testid="output_38" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="1C3D3C91" data-testid="output_39" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5668 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="B8847B1A" data-testid="output_40" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D87FF9AB" data-testid="output_41" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5674 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="B822A190" data-testid="output_42" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="356B817C" data-testid="output_43" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5680 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="39C11B5E" data-testid="output_44" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="756CF658" data-testid="output_45" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5686 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="661F69D1" data-testid="output_46" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="F9CE4E2C" data-testid="output_47" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5691 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="10FFA469" data-testid="output_48" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="6DC81488" data-testid="output_49" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5697 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="1E1DB440" data-testid="output_50" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="5B7A988D" data-testid="output_51" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5703 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="4143B426" data-testid="output_52" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="4F012A68" data-testid="output_53" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5708 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="E4CD9BB9" data-testid="output_54" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="067B4422" data-testid="output_55" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5714 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="C3717DA6" data-testid="output_56" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D2FFA638" data-testid="output_57" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Growth rate 0.5720 :</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D5A9AC24" data-testid="output_58" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  No	Perc	      Name	      Min	      Max</div></div></div></div></div><div  class = 'S6'><span>Plot the results to visualize the difference (see also Figure 2 in ref. [1]):</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >grComV = result.GRmax * options.optGRpercent / 100;  </span><span style="color: rgb(2, 128, 9);">% vector of growth rates tested</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >lgLabel = {</span><span style="color: rgb(170, 4, 249);">'{\itEc1 }'</span><span >;</span><span style="color: rgb(170, 4, 249);">'{\itEc2 }'</span><span >;</span><span style="color: rgb(170, 4, 249);">'{\itEc3 }'</span><span >;</span><span style="color: rgb(170, 4, 249);">'{\itEc4 }'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >col = [235 135 255; 0 235 0; 255 0 0; 95 135 255 ]/255;  </span><span style="color: rgb(2, 128, 9);">% color</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >f = figure;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% SteadyCom</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >subplot(2, 1, 1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">on</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >x = [grComV(:); flipud(grComV(:))];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >j = 1:4</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    y = [fvaComMin(j, :), fliplr(fvaComMax(j, :))];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    p(j, 1) = plot(x(~isnan(y)), y(~isnan(y)), </span><span style="color: rgb(170, 4, 249);">'LineWidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    p(j, 1).Color = col(j, :);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >tl(1) = title(</span><span style="color: rgb(170, 4, 249);">'\underline{SteadyCom}'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Interpreter'</span><span >, </span><span style="color: rgb(170, 4, 249);">'latex'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >tl(1).Position = [0.7 1.01 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ax(1) = gca;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ax(1).XTick = 0.66:0.02:0.74;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ax(1).YTick = 0:0.2:1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >xlim([0.66 0.74])</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ylim([0 1])</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >lg = legend(lgLabel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >lg.Box = </span><span style="color: rgb(170, 4, 249);">'off'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >yl(1) = ylabel(</span><span style="color: rgb(170, 4, 249);">'Relative abundance'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >xl(1) = xlabel(</span><span style="color: rgb(170, 4, 249);">'Community growth rate (h^{-1})'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% FBA</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >grFBAV = grFBA * optGRpercentFBA / 100;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >x = [grFBAV(:); flipud(grFBAV(:))];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >subplot(2, 1, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">on</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% plot j=1:2 only because 3:4 overlap with 1:2</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >j = 1:2</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    y = [fvaFBAMin(j, :), fliplr(fvaFBAMax(j, :))] ./ x';</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(2, 128, 9);">% it is possible some values &gt; 1 because the total biomass produced is</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(2, 128, 9);">% only bounded below when calling fluxVariability. Would be strictly</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(2, 128, 9);">% equal to 1 if sum(biomass) = optGRpercentFBA(jGr) * grFBA is constrained. Treat them as 1.</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    y(y&gt;1) = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    p(j, 2)= plot(x(~isnan(y)), y(~isnan(y)), </span><span style="color: rgb(170, 4, 249);">'LineWidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    p(j, 2).Color = col(j, :);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >tl(2) = title(</span><span style="color: rgb(170, 4, 249);">'\underline{Joint FBA}'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Interpreter'</span><span >, </span><span style="color: rgb(170, 4, 249);">'latex'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >tl(2).Position = [0.55 1.01 0];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ax(2) = gca;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ax(2).XTick = 0.52:0.02:0.58;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ax(2).YTick = 0:0.2:1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >xlim([0.52 0.58])</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ylim([0 1])</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >xl(2) = xlabel(</span><span style="color: rgb(170, 4, 249);">'Community growth rate (h^{-1})'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >yl(2) = ylabel(</span><span style="color: rgb(170, 4, 249);">'Relative abundance'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ax(1).Position = [0.1 0.6 0.5 0.32];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ax(2).Position = [0.1 0.1 0.5 0.32];</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >lg.Position = [0.65 0.65 0.1 0.27];</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="91A9B4BC" data-testid="output_59" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S6'><span>The direct use of FVA compared to FVA under the SteadyCom framework gives very little information on the organism's abundance. The ranges for almost all growth rates span from 0 to 1. In contrast, </span><span style=' font-family: monospace;'>SteadyComFVA</span><span> returns results with the expected co-existence of all four mutants. When the growth rates get closer to the maximum, the ranges shrink to unique values.</span></div><div  class = 'S2'><span></span></div><h2  class = 'S7'><span>Analyze Pairwise Relationship Using SteadyComPOA</span></h2><div  class = 'S2'><span>Now we would like to see at a given growth rate, how the abundance of an organism influences the abundance of another organism. We check this by iteratively fixing the abundance of an organism at a level (independent variable) and optimizing for the maximum and minimum allowable abundance of another organism (dependent variable). This is what </span><span style=' font-family: monospace;'>SteadyComPOA</span><span> does.</span></div><div  class = 'S2'><span>Set up the option structure and call </span><span style=' font-family: monospace;'>SteadyComPOA</span><span>. </span><span style=' font-family: monospace;'>Nstep</span><span> is an important parameter to designate how many intermediate steps are used or which values between the min and max values of the independent variable are used for optimizing the dependent variable. </span><span style=' font-family: monospace;'>savePOA</span><span> options must be supplied with a non-empty string or a default name will be used for saving the POA results. By default, the function analyzes all possible pairs in </span><span style=' font-family: monospace;'>options.rxnNameList</span><span>. To analyze only particular pairs, use </span><span style=' font-family: monospace;'>options.pairList</span><span>. See </span><span style=' font-family: monospace;'>help SteadyComPOA </span><span>for more details.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >options.savePOA = [</span><span style="color: rgb(170, 4, 249);">'POA' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'EcCom'</span><span >];  </span><span style="color: rgb(2, 128, 9);">% directory and fila name for saving POA results</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.optGRpercent = [99 90 70 50];  </span><span style="color: rgb(2, 128, 9);">% analyze at these percentages of max. growth rate</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Nstep is the number of intermediate steps that the independent variable will take different values</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% or directly the vector of values, e.g. Nsetp = [0, 0.5, 1] implies fixing the independent variable at the minimum,</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% 50% from the min to the max and the maximum value respectively to find the attainable range of the dependent variable.</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Here use small step sizes when getting close to either ends of the flux range</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >a = 0.001*(1000.^((0:14)/14));</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.Nstep = sort([a (1-a)]);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >[POAtable, fluxRange, Stat, GRvector] = SteadyComPOA(EcCom, options);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="F14AFBD9" data-testid="output_60" data-width="420" data-height="59" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Already finished. Results were already saved to POA/EcCom_GR0.73.mat
Already finished. Results were already saved to POA/EcCom_GR0.66.mat
Already finished. Results were already saved to POA/EcCom_GR0.52.mat
Already finished. Results were already saved to POA/EcCom_GR0.37.mat</div></div></div></div></div><div  class = 'S6'><span>POAtable is a </span><span style=' font-style: italic;'>n</span><span>-by-</span><span style=' font-style: italic;'>n</span><span> cell if there are </span><span style=' font-style: italic;'>n</span><span> targets in </span><span style=' font-family: monospace;'>options.rxnNameList</span><span>. </span><span style=' font-family: monospace;'>POAtable{i, i}</span><span> is a </span><span style=' font-style: italic;'>Nstep</span><span>-by-1-by-</span><span style=' font-style: italic;'>Ngr</span><span> matrix where </span><span style=' font-style: italic;'>Nstep </span><span>is the number of intermediate steps detemined by </span><span style=' font-family: monospace;'>options.Nstep</span><span> and </span><span style=' font-style: italic;'>Ngr </span><span>is the number of growth rates analyzed. </span><span style=' font-family: monospace;'>POAtable{i, i}(:, :, k)</span><span> is the values at which the </span><span style=' font-style: italic;'>i</span><span>-th target is fixed for the community growing at the growth rate </span><span style=' font-family: monospace;'>GRvector(k)</span><span>. POAtable{i, j} is a </span><span style=' font-style: italic;'>Nstep</span><span>-by-2-by-</span><span style=' font-style: italic;'>Ngr</span><span> matrix where </span><span style=' font-family: monospace;'>POAtable{i, j}(:, 1, k)</span><span> and </span><span style=' font-family: monospace;'>POAtable{i, j}(:, 2, k)</span><span> are respectively the min. and max. values of the </span><span style=' font-style: italic;'>j</span><span>-th target when fixing the </span><span style=' font-style: italic;'>i</span><span>-th target at the corresponding values in </span><span style=' font-family: monospace;'>POAtable{i, i}(:, :, k)</span><span>. </span><span style=' font-family: monospace;'>fluxRange </span><span>contains the min. and max. values for each target (found by calling </span><span style=' font-family: monospace;'>SteadyComFVA</span><span>). </span><span style=' font-family: monospace;'>Stat </span><span>is a </span><span style=' font-style: italic;'>n</span><span>-by-</span><span style=' font-style: italic;'>n-by-Ngr</span><span> structure array, each containing two fields: </span><span style=' font-family: monospace;'>*.cor</span><span>, the correlatiion coefficient between the max/min values of the dependent variable and the independent variable, and </span><span style=' font-family: monospace;'>*.r2</span><span>, the R-squred of linear regression. They are also outputed in the command window during computation. All the computed results are also saved in the folder 'POA' starting with the name 'EcCom', followed by 'GRxxxx' denoting the growth rate at which the analysis is performed.</span></div><div  class = 'S2'><span>Plot the results (see also Figure 3 in ref. [1]):</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >nSp = 4;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >spLab = {</span><span style="color: rgb(170, 4, 249);">'{\it Ec1 }'</span><span >;</span><span style="color: rgb(170, 4, 249);">'{\it Ec2 }'</span><span >;</span><span style="color: rgb(170, 4, 249);">'{\it Ec3 }'</span><span >;</span><span style="color: rgb(170, 4, 249);">'{\it Ec4 }'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >mark = {</span><span style="color: rgb(170, 4, 249);">'A'</span><span >, </span><span style="color: rgb(170, 4, 249);">'B'</span><span >, </span><span style="color: rgb(170, 4, 249);">'D'</span><span >, </span><span style="color: rgb(170, 4, 249);">'C'</span><span >, </span><span style="color: rgb(170, 4, 249);">'E'</span><span >, </span><span style="color: rgb(170, 4, 249);">'F'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >nPlot = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >j = 1:nSp</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >k = 1:nSp</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >k &gt; j</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            nPlot = nPlot + 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k) = subplot(nSp-1, nSp-1, (k - 2) * (nSp - 1) + j);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            hold </span><span style="color: rgb(170, 4, 249);">on</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            </span><span style="color: rgb(14, 0, 255);">for </span><span >p = 1:size(POAtable{1, 1}, 3)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >                x = [POAtable{j, j}(:, :, p);POAtable{j, j}(end:-1:1, :, p);</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >                    POAtable{j, j}(1, 1, p)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >                y = [POAtable{j, k}(:, 1, p);POAtable{j, k}(end:-1:1, 2, p);</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >                        POAtable{j, k}(1, 1, p)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >                plot(x(~isnan(y)), y(~isnan(y)), </span><span style="color: rgb(170, 4, 249);">'LineWidth'</span><span >, 2)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            xlim([0.001 1])</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ylim([0.001 1])</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).XScale = </span><span style="color: rgb(170, 4, 249);">'log'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).YScale = </span><span style="color: rgb(170, 4, 249);">'log'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).XTick = [0.001 0.01 0.1 1];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).YTick = [0.001 0.01 0.1 1];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).YAxis.MinorTickValues=[];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).XAxis.MinorTickValues=[];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).TickLength = [0.03 0.01];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            xlabel(spLab{j});</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ylabel(spLab{k});</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            tx(j, k) = text(10^(-5), 10^(0.1), mark{nPlot}, </span><span style="color: rgb(170, 4, 249);">'FontSize'</span><span >, 12, </span><span style="color: rgb(170, 4, 249);">'FontWeight'</span><span >, </span><span style="color: rgb(170, 4, 249);">'bold'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >lg = legend(strcat(strtrim(cellstr(num2str(options.optGRpercent(:)))), </span><span style="color: rgb(170, 4, 249);">'%'</span><span >));</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >lg.Position = [0.7246 0.6380 0.1700 0.2015];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >lg.Box=</span><span style="color: rgb(170, 4, 249);">'off'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >subplot(3, 3, 3, </span><span style="color: rgb(170, 4, 249);">'visible'</span><span >, </span><span style="color: rgb(170, 4, 249);">'off'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >t = text(0.2, 0.8, {</span><span style="color: rgb(170, 4, 249);">'% maximum'</span><span >;</span><span style="color: rgb(170, 4, 249);">'growth rate'</span><span >});</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >j = 1:nSp</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >k = 1:nSp</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >k&gt;j</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).Position = [0.15 + (j - 1) * 0.3, 0.8 - (k - 2) * 0.3, 0.16, 0.17];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >            ax(j, k).Color = </span><span style="color: rgb(170, 4, 249);">'none'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >        </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="6BF5B4B8" data-testid="output_61" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S6'><span>There are two patterns observed. The two pairs showing negative correlations, namely Ec1 vs Ec4 (panel D) and Ec2 vs Ec3 (panel C) are indeed competing for the same amino acids with each other (Ec1 and Ec4 competing for Lys and Met; Ec2 and Ec4 competing for Arg and Phe). Each of the other pairs showing positive correlations are indeed the cross feeding pairs, e.g., Ec1 and Ec2 (panel A) cross feeding on Arg and Lys. See ref. [1] for more detailed discussion.</span></div><div  class = 'S2'><span></span></div><h2  class = 'S7'><span>Parallelization and Timing</span></h2><div  class = 'S2'><span>SteadyCom in general can be finished within 20 iterations, i.e. solving 20 LPs (usually faster if using Matlab </span><span style=' font-family: monospace;'>fzero</span><span>) for an accuracy of 1e-6 for the maximum community growth rate. The actual computation time depends on the size of the community metabolic network. The current </span><span style=' font-family: monospace;'>EcCom</span><span> model has 6971 metabolites and 9831 reactions. It took 18 seconds for a MacBook Pro with 2.5 GHz Intel Core i5, 4 GB memory running Matlab R2016b and Cplex 12.7.1.</span></div><div  class = 'S2'><span>Since the FVA and POA analysis can be time-consuming for large models with a large number of reactions to be analyzed, SteadyComFVA and SteadyComPOA support parrallelization using the Matlab Distributed Computing Toolbox (</span><span style=' font-family: monospace;'>parfor</span><span> for SteadyComFVA and </span><span style=' font-family: monospace;'>spmd</span><span> for SteadyComPOA). </span></div><div  class = 'S2'><span>Test SteadyComFVA with 2 threads:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >options.rxnNameList = EcCom.rxns(1:100);  </span><span style="color: rgb(2, 128, 9);">% test FVA for the first 50 reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.optGRpercent = 99;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.algorithm = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.threads = 1;  </span><span style="color: rgb(2, 128, 9);">% test single-thread computation first</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.verbFlag = 0;  </span><span style="color: rgb(2, 128, 9);">% no verbose output</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >tic;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >[minF1, maxF1] = SteadyComFVA(EcCom, options);  </span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >t1 = toc;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">if </span><span >isempty(gcp(</span><span style="color: rgb(170, 4, 249);">'nocreate'</span><span >))</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    parpool(2);  </span><span style="color: rgb(2, 128, 9);">% start a parallel pool</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="FDA6DDC8" data-testid="output_62" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Starting parallel pool (parpool) using the 'local' profile ... connected to 2 workers.</div></div></div></div><div class="inlineWrapper"><div  class = 'S13'><span style="white-space: pre"><span >options.threads = 2;  </span><span style="color: rgb(2, 128, 9);">% two threads (0 to use all available workers)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >tic;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >[minF2, maxF2] = SteadyComFVA(EcCom, options);  </span><span style="color: rgb(2, 128, 9);">% test single-thread computation first</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >t2 = toc;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'Maximum difference between the two solutions: %.4e\n'</span><span >, max(max(abs(minF1 - minF2)), max(abs(maxF1 - maxF2))));</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="B570A227" data-testid="output_63" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Maximum difference between the two solutions: 9.9257e-09</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'\nSingle-thread computation: %.0f sec\nTwo-thread computation: %.0f sec\n'</span><span >, t1, t2);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="02BF795B" data-testid="output_64" data-width="420" data-height="31" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Single-thread computation: 96 sec
Two-thread computation: 91 sec</div></div></div></div></div><div  class = 'S6'><span>If there are many reactions to be analyzed, use </span><span style=' font-family: monospace;'>options.saveFVA</span><span> to give a relative path for saving the intermediate results. Even though the computation is interrupted, by calling </span><span style=' font-family: monospace;'>SteadyComFVA</span><span> with the same </span><span style=' font-family: monospace;'>options.saveFVA</span><span>, the program will detect previously saved results and continued from there.</span></div><div  class = 'S2'><span>Test SteadyComPOA with 2 threads:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >options.rxnNameList = EcCom.rxns(find(abs(result.flux) &gt; 1e-2, 6));</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.savePOA = </span><span style="color: rgb(170, 4, 249);">'POA/EcComParallel'</span><span >;  </span><span style="color: rgb(2, 128, 9);">% save with a new name</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.verbFlag = 3;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.threads = 2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.Nstep = 5;  </span><span style="color: rgb(2, 128, 9);">% use a smaller number of steps for test</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >tic;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >[POAtable1, fluxRange1] = SteadyComPOA(EcCom, options);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="804FD91F" data-testid="output_65" data-width="420" data-height="885" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Find maximum community growth rate..
Model feasible at maintenance. Time elapsed: 1 / 1 sec
Iter        LB   To test        UB  Time elapsed (iteration/total)
   1  0.000000  0.500000       Inf  0 / 1 sec
   2  0.500000  0.721279       Inf  6 / 7 sec
   3  0.721279  0.735372       Inf  0 / 7 sec
   4  0.735372  0.742726       Inf  0 / 8 sec
 
 Func-count    x          f(x)             Procedure
    2        0.735372  -0.000807615        initial
    3        0.735378   -0.00079987        interpolation
    4         0.73599  -1.26127e-06        interpolation
    5         0.73599  -1.26127e-06        interpolation
 
Zero found in the interval [0.735372, 0.742726]
Maximum community growth rate: 0.735990 (abs. error &lt; 1e-06).	Time elapsed: 26 sec

FVA for 6 sets of fluxes/biomass at growth rate 0.728630 :

Thread 1:	33.33% finished. 2017-07-21 13:56:18
Thread 2:	33.33% finished. 2017-07-21 13:56:18
Thread 1:	66.67% finished. 2017-07-21 13:56:20
Thread 2:	66.67% finished. 2017-07-21 13:56:20
Thread 1:	100.00% finished. 2017-07-21 13:56:21
Thread 2:	100.00% finished. 2017-07-21 13:56:21

POA for 15 pairs of reactions at growth rate 0.728630
Start from #1 Ec13HAD100 vs #2 Ec13HAD120.
           Rxn1           Rxn2    corMin        r2    corMax        r2   Time
POA in parallel...
Lab 2: 
       Ec13HAD120     Ec13HAD160    0.0956    0.5000   -0.8431    0.9667   2017-07-21 13:57:45
Lab 1: 
       Ec13HAD100     Ec13HAD120    0.5755    0.3373    0.7927    0.4005   2017-07-21 13:58:23
Lab 2: 
       Ec13HAD121     Ec13HAD140   -0.0837    0.5000   -0.3890    0.0784   2017-07-21 13:59:32
Lab 1: 
       Ec13HAD100     Ec13HAD121    0.2429    0.7227    0.4245    0.2168   2017-07-21 14:00:44
Lab 2: 
       Ec13HAD121     Ec13HAD141    0.9997    1.0000    1.0000    1.0000   2017-07-21 14:01:18
Lab 1: 
       Ec13HAD100     Ec13HAD140   -0.0915    0.4667   -0.1144    1.0000   2017-07-21 14:01:54
Lab 2: 
       Ec13HAD121     Ec13HAD160   -0.0837    0.5000   -0.2478    0.0302   2017-07-21 14:02:33
       Ec13HAD140     Ec13HAD141   -0.0197    0.1369   -0.6518    0.9578   2017-07-21 14:04:17
Lab 1: 
       Ec13HAD100     Ec13HAD141    0.2429    0.7226    0.4245    0.2447   2017-07-21 14:04:52
       Ec13HAD100     Ec13HAD160    0.0000       NaN    1.8482    0.4493   2017-07-21 14:05:44
       Ec13HAD120     Ec13HAD121   -0.0922    0.3440   -0.5288    0.9995   2017-07-21 14:07:32
Lab 2: 
       Ec13HAD140     Ec13HAD160    0.1842    0.8929   -1.0433    0.9735   2017-07-21 14:08:04
       Ec13HAD141     Ec13HAD160   -0.0837    0.5000   -0.2478    0.0302   2017-07-21 14:09:11
Lab 1: 
       Ec13HAD120     Ec13HAD140   -0.0000       NaN   -1.4156    1.0000   2017-07-21 14:09:25
Lab 2: 
  Current loop finished. Stop other workers...
  All workers have ceased. Redistributing...
Lab 1: 
       Ec13HAD120     Ec13HAD141   -0.0402    0.1302   -0.6122    0.9816   2017-07-21 14:10:29
Lab 2: 
  Current loop finished. Stop other workers...
  All workers have ceased. Redistributing...
Finished. Save final results to POA/EcComParallel_GR0.73.mat</div></div></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >t3 = toc;</span></span></div></div></div><div  class = 'S6'><span>The parallelization code uses </span><span style=' font-family: monospace;'>spmd</span><span> and will redistribute jobs once any of the workers has finished to maximize the computational efficiency.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >options.savePOA = </span><span style="color: rgb(170, 4, 249);">'POA/EcComSingeThread'</span><span >;  </span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.threads = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >tic;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >[POAtable2, fluxRange2] = SteadyComPOA(EcCom, options);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="8CB0E405" data-testid="output_66" data-width="420" data-height="633" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Find maximum community growth rate..
Model feasible at maintenance. Time elapsed: 1 / 1 sec
Iter        LB   To test        UB  Time elapsed (iteration/total)
   1  0.000000  0.500000       Inf  0 / 1 sec
   2  0.500000  0.721279       Inf  5 / 6 sec
   3  0.721279  0.735372       Inf  0 / 6 sec
   4  0.735372  0.742726       Inf  0 / 6 sec
 
 Func-count    x          f(x)             Procedure
    2        0.735372  -0.000807615        initial
    3        0.735378   -0.00079987        interpolation
    4         0.73599  -1.26127e-06        interpolation
    5         0.73599  -1.26127e-06        interpolation
 
Zero found in the interval [0.735372, 0.742726]
Maximum community growth rate: 0.735990 (abs. error &lt; 1e-06).	Time elapsed: 24 sec

FVA for 6 sets of fluxes/biomass at growth rate 0.728630 :
  No	   %	      Name	      Min	      Max
   1	  17	Ec13HAD100	 0.052591	 0.217439
   2	  33	Ec13HAD120	 0.000000	 0.262936
   3	  50	Ec13HAD121	 0.022231	 0.202541
   4	  67	Ec13HAD140	 0.000000	 0.243774
   5	  83	Ec13HAD141	 0.022231	 0.202541
   6	 100	Ec13HAD160	 0.000000	 0.251518

POA for 15 pairs of reactions at growth rate 0.728630
Start from #1 Ec13HAD100 vs #2 Ec13HAD120.
           Rxn1           Rxn2    corMin        r2    corMax        r2   Time
     Ec13HAD100     Ec13HAD120    0.5755    0.3373    0.7927    0.4005   2017-07-21 14:11:54
     Ec13HAD100     Ec13HAD121    0.2429    0.7227    0.4245    0.2168   2017-07-21 14:13:16
     Ec13HAD100     Ec13HAD140   -0.0915    0.4667   -0.1144    1.0000   2017-07-21 14:13:54
     Ec13HAD100     Ec13HAD141    0.2429    0.7226    0.4245    0.2447   2017-07-21 14:15:10
     Ec13HAD100     Ec13HAD160    0.0000       NaN    1.8482    0.4493   2017-07-21 14:15:39
     Ec13HAD120     Ec13HAD121   -0.0922    0.3440   -0.5288    0.9995   2017-07-21 14:16:30
     Ec13HAD120     Ec13HAD140   -0.0000       NaN   -1.4156    1.0000   2017-07-21 14:17:36
     Ec13HAD120     Ec13HAD141    0.0637    1.0000   -0.6611    0.9793   2017-07-21 14:18:38
     Ec13HAD120     Ec13HAD160    0.1435    0.6000   -0.8448    0.9673   2017-07-21 14:18:48
     Ec13HAD121     Ec13HAD140   -0.0837    0.5000   -0.3890    0.0784   2017-07-21 14:19:34
     Ec13HAD121     Ec13HAD141    0.9997    1.0000    1.0000    1.0000   2017-07-21 14:20:18
     Ec13HAD121     Ec13HAD160   -0.0837    0.5000   -0.2478    0.0302   2017-07-21 14:20:44
     Ec13HAD140     Ec13HAD141   -0.0026    0.0014   -0.6518    0.9589   2017-07-21 14:21:16
     Ec13HAD140     Ec13HAD160    0.1547    0.6000   -0.9028    1.0000   2017-07-21 14:22:06
     Ec13HAD141     Ec13HAD160   -0.0837    0.4667   -0.2437    0.0293   2017-07-21 14:22:51
Finished. Save final results to POA/EcComSingeThread_GR0.73.mat</div></div></div></div><div class="inlineWrapper"><div  class = 'S13'><span style="white-space: pre"><span >t4 = toc;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >dev = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >i = 1:size(POAtable1, 1)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >j = i:size(POAtable1, 2)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >        dev = max(max(max(abs(POAtable1{i, j} - POAtable2{i, j}))));</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >        dev = max(dev, max(max(abs(fluxRange1 - fluxRange2))));</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'Maximum difference between the two solutions: %.4e\n'</span><span >, dev);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="E70E04DF" data-testid="output_67" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Maximum difference between the two solutions: 1.7043e-09</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'\nSingle-thread computation: %.0f sec\nTwo-thread computation: %.0f sec\n'</span><span >, t4, t3);</span></span></div><div  class = 'S5'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="B51A5CF2" data-testid="output_68" data-width="420" data-height="31" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Single-thread computation: 742 sec
Two-thread computation: 879 sec</div></div></div></div></div><div  class = 'S6'><span>The advantage will be more significant for more targets to analyzed and more threads used. Similar to </span><span style=' font-family: monospace;'>SteadyComFVA</span><span>, </span><span style=' font-family: monospace;'>SteadyComPOA</span><span> also supports continuation from previously interrupted computation by calling with the same </span><span style=' font-family: monospace;'>options.savePOA</span><span>.</span></div><div  class = 'S2'><span> </span></div><h2  class = 'S7'><span>REFERENCES</span></h2><div  class = 'S2'><span style=' font-style: italic;'>[1] </span><span>Chan SHJ, Simons MN, Maranas CD (2017) SteadyCom: Predicting microbial abundances while ensuring community stability. PLoS Comput Biol 13(5): e1005539. https://doi.org/10.1371/journal.pcbi.1005539</span></div><div  class = 'S2'><span style=' font-style: italic;'>[2] </span><span>Khandelwal RA, Olivier BG, Röling WFM, Teusink B, Bruggeman FJ (2013) Community Flux Balance Analysis for Microbial Consortia at Balanced Growth. PLoS ONE 8(5): e64567. https://doi.org/10.1371/journal.pone.0064567</span></div><div  class = 'S2'></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Analyze Steady-State Community COBRA Models
%% Author(s): Siu Hung Joshua Chan, Department of Chemical Engineering, The Pennsylvania State University
%% Reviewer(s): Almut Heinken, Luxembourg Centre for Systems Biomedicine, University of Luxembourg
% 
%% INTRODUCTION
% This tutorial demonstrates the use of SteadyCom to analyze a multi-organism 
% COBRA model (e.g., for a microbial community) at a community steady-state [1]. 
% Compared to the direct extension of flux balance analysis (FBA) which simply 
% treats a community model as a multi-compartment model, SteadyCom explicitly 
% introduces the biomass variables to describe the relationships between biomass, 
% biomass production rate, growth rate and fluxes. SteadyCom also assumes the 
% existence of a time-averaged population steady-state for a stable microbial 
% community which in turn implies a time-averaged constant growth rate across 
% all members. SteadyCom is equivalent to the reformulation of the earlier community 
% flux balance analysis (cFBA) [2] with significant computational advantage. SteadyCom 
% computes the maximum community growth rate by solving the follow optimization 
% problem:
% 
% $$\begin{array}{ll}\max & \ \mu\\ \\\text{s.t.} & \sum\limits_{j\in \textbf{J}^k}S^k_{ij}V^k_j=0, 
% & \forall i \in \textbf{I}^k, k \in \textbf{K}\\ & LB^k_jX^k\leq V^k_j\leq UB^k_jX^k, 
% & \forall j \in \textbf{J}^k, k \in \textbf{K} \\& \sum\limits_{k \in \textbf{K}}V^k_{ex(i)} 
% + u^{com}_i\geq 0, & \forall i \in \textbf{I}^{com} \\& V^k_{biomass} = X^k\mu, 
% & \forall k \in \textbf{K} \\& \sum\limits_{k \in \textbf{K}}X^k = 1 \\& X^k,\quad 
% \mu \geq 0, & \forall k \in \textbf{K} \\& V^k_j \in \Re, & \forall j \in \textbf{J}^k, 
% k \in \textbf{K} \end{array}$$
% 
% where $S^k_{ij}$ is the stoichiometry of metabolite _i_ in reaction _j_ for 
% organism _k_, $V^k_j$, $LB^k_j$ and $UB^k_j$ are respectively the flux (in mmol/h), 
% lower bound (in mmol/h/gdw) and upper bound (in mmol/h/gdw) for reaction _j_ 
% for organism _k_, $u^{com}_i$ is the community uptake bound for metabolite _i_, 
% $X^k$ is the biomass (in gdw) of organism _k_, $\mu$ is the community growth 
% rate, $\textbf{I}^k$ is the set of metabolites of organism _k_, $\textbf{I}^{com}$ 
% is the set of community metabolites in the community exchange space, $\textbf{J}^k$ 
% is the set of reactions for organism k, $\textbf{K}$ is the set of organisms 
% in the community, and $ex(i) \in \textbf{J}^k$ is the exchange reaction in organism 
% _k_ for extracellular metabolite _i_. See ref. [1] for the derivation and detailed 
% explanation.
% 
% Throughout the tutorial, using a hypothetical model of four _E. coli_ mutants 
% auxotrophic for amino acids, we will demonstrate the three different functionalities 
% of the module: (1) computing the maximum community growth rate using the function 
% SteadyCom.m, (2) performing flux variability analysis under a given community 
% growth rate using SteadyComFVA.m, and (3) analyzing the pairwise relationship 
% between flux/biomass variables using a technique similar to Pareto-optimal analysis 
% by calling the function SteadyComPOA.m
%% EQUIPMENT SETUP
% If necessary, initialise the cobra toolbox and select a solver by running:

initCobraToolbox(false) % false, as we don't want to update
%% 
% All SteadyCom functions involve only solving linear programming problems. 
% Any solvers supported by the COBRA toolbox will work. But SteadyCom contains 
% specialized codes for IBM ILOG Cplex which was tested to run significantly faster 
% for SteadyComFVA and SteadyComPOA for larger problems through calling the Cplex 
% object in Matlab directly. 
% 
% Please note that parallelization requires a working installation of the Parallel 
% Computing Toolbox.

changeCobraSolver('ibm_cplex', 'LP');
%% PROCEDURE
%% Model Construction
% Load the _E. coli_ iAF1260 model in the COBRA toolbox.

global CBTDIR
iAF1260 = readCbModel([CBTDIR filesep 'test' filesep 'models' filesep 'iAF1260.mat']);
%% 
% Polish the model a little bit:

% convert the compartment format from e.g., '_c' to '[c]'
iAF1260.mets = regexprep(iAF1260.mets, '_([^_]+)$', '\[$1\]');
% make all empty cells in cell arrays to be empty string
fieldToBeCellStr = {'metFormulas'; 'genes'; 'grRules'; 'metNames'; 'rxnNames'; 'subSystems'};
for j = 1:numel(fieldToBeCellStr)
    iAF1260.(fieldToBeCellStr{j})(cellfun(@isempty, iAF1260.(fieldToBeCellStr{j}))) = {''};
end
%% 
% Add a methionine export reaction to allow the export of methionine.

iAF1260 = addReaction(iAF1260,{'METt3pp',''},'met__L[c] + h[c] => met__L[p] + h[p]');
%% 
% Reactions essential for amino acid autotrophy:

argH = {'ARGSL'};  % essential for arginine biosynthesis
lysA = {'DAPDC'};  % essential for lysine biosynthesis
metA = {'HSST'};  % essential for methionine biosynthesis
ilvE = {'PPNDH'};  % essential for phenylalanine biosynthesis
%% 
% Reactions essential for exporting amino acids:

argO = {'ARGt3pp'};  % Evidence for an arginine exporter encoded by yggA (argO) that is regulated by the LysR-type transcriptional regulator ArgP in Escherichia coli.
lysO = {'LYSt3pp'};  % Distinct paths for basic amino acid export in Escherichia coli: YbjE (LysO) mediates export of L-lysine
yjeH = {'METt3pp'};  % YjeH is a novel L-methionine and branched chain amino acids exporter in Escherichia coli
yddG = {'PHEt2rpp'};  % YddG from Escherichia coli promotes export of aromatic amino acids.
%% 
% Now make four copies of the model with auxotrophy for different amino acids 
% and inability to export amino acids:

% auxotrophic for Lys and Met, not exporting Phe
Ec1 = iAF1260;
Ec1 = changeRxnBounds(Ec1, [lysA; metA; yddG], 0, 'b');
% auxotrophic for Arg and Phe, not exporting Met
Ec2 = iAF1260;
Ec2 = changeRxnBounds(Ec2, [argH; yjeH; ilvE], 0, 'b');
% Auxotrophic for Arg and Phe, not exporting Lys
Ec3 = iAF1260;
Ec3 = changeRxnBounds(Ec3, [argH; lysO; ilvE], 0, 'b');
% Auxotrophic for Lys and Met, not exporting Arg
Ec4 = iAF1260;
Ec4 = changeRxnBounds(Ec4, [argO; lysA; metA], 0, 'b');
%% 
% Now none of the four organisms can grow alone and they must cross feed each 
% other to survive. See Figure 1 in ref. <http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005539 
% [1]> for the visualization of the community. 
% 
% Get the extracellular metabolites, the corresponding exchange reactions and 
% the uptake rates for the _E. coli_ model, which are used later to constrain 
% the community model:

% extracellular metabolites (met[e])
metEx = strcmp(getCompartment(iAF1260.mets),'e');
% the corresponding exchange reactions
rxnExAll = find(sum(iAF1260.S ~= 0, 1) == 1);
[rxnEx, ~] = find(iAF1260.S(metEx, rxnExAll)');  % need to be in the same order as metEx
rxnEx = rxnExAll(rxnEx);
% exchange rate
lbEx = iAF1260.lb(rxnEx);
%% 
% Create a community model with the four _E. coli_ tagged as 'Ec1', 'Ec2', 'Ec3', 
% 'Ec4' respectively by calling |createMultipleSpeciesModel|. 

nameTagsModel = {'Ec1'; 'Ec2'; 'Ec3'; 'Ec4'};
EcCom = createMultipleSpeciesModel({Ec1; Ec2; Ec3; Ec4}, nameTagsModel);
EcCom.csense = char('E' * ones(1,numel(EcCom.mets)));  % correct the csense
clear Ec1 Ec2 Ec3 Ec4
%% 
% The model |EcCom| contains a community compartment denoted by |[u]| to allow 
% exchange between organisms. Each organism-specific reaction/metabolite is prepended 
% with the corresponding tag.
% 
% Retreive the names and ids for organism/community exchange reactions/metabolites 
% which are necessary for computation:

[EcCom.infoCom, EcCom.indCom] = getMultiSpeciesModelId(EcCom, nameTagsModel);
disp(EcCom.infoCom);
%% 
% |EcCom.infoCom| contains reaction/metabolite names (from |EcCom.rxns|/|EcCom.mets|) 
% for the community exchange reactions (|*.EXcom|), organism-community exchange 
% reactions (|*.EXsp|), community metabolites (|*.Mcom|), organism-specific extracellular 
% metabolite (|*.Msp|). If a host model is specified, there will also be non-empty 
% |*.EXhost| and |*.Mhost| for the host-specific exchange reactions and metabolites. 
% The fields |*.rxnSps|/|*.metSps| give information on which organism a reaction/metabolite 
% belongs to.
% 
% |indCom| has the same structure as |infoCom| but contains the indices rather 
% than names. |infoCom| and |indCom| are attached as fields of the model |EcCom| 
% because SteadyCom requires this information from the input model for computation. 
% Incorporate also the names and indices for the biomass reactions which are necessary 
% for computing growth:

rxnBiomass = strcat(nameTagsModel, 'BIOMASS_Ec_iAF1260_core_59p81M');  % biomass reaction names
rxnBiomassId = findRxnIDs(EcCom, rxnBiomass);  % ids
EcCom.infoCom.spBm = rxnBiomass;  % .spBm for organism biomass reactions
EcCom.indCom.spBm = rxnBiomassId;
%% 
% 
%% Finding Maximum Growth Rate Using SteadyCom
% Set community and organism-specific uptake rates to be the same as in the 
% orginal iAF1260 model:

[yn, id] = ismember(strrep(iAF1260.mets(metEx), '[e]', '[u]'), EcCom.infoCom.Mcom);  % map the metabolite name
assert(all(yn));  % must be a 1-to-1 mapping
EcCom.lb(EcCom.indCom.EXcom(:,1)) = lbEx(id);  % assign community uptake bounds
EcCom.ub(EcCom.indCom.EXcom(:,1)) = 1e5;
EcCom.lb(EcCom.indCom.EXsp) = repmat(lbEx(id), 1, 4);  % assign organism-specific uptake bounds
%% 
% Set maximum allowed organism-specific uptake rates for the cross-feeding amino 
% acids:

% only allow to take up the amino acids that one is auxotrophic for
exRate = 1;  % maximum uptake rate for cross feeding AAs
% Ec1
EcCom = changeRxnBounds(EcCom, {'Ec1IEX_arg__L[u]tr'; 'Ec1IEX_phe__L[u]tr'}, 0, 'l');
EcCom = changeRxnBounds(EcCom, {'Ec1IEX_met__L[u]tr'; 'Ec1IEX_lys__L[u]tr'}, -exRate, 'l');
% Ec2
EcCom = changeRxnBounds(EcCom, {'Ec2IEX_arg__L[u]tr'; 'Ec2IEX_phe__L[u]tr'}, -exRate, 'l');
EcCom = changeRxnBounds(EcCom, {'Ec2IEX_met__L[u]tr'; 'Ec2IEX_lys__L[u]tr'}, 0, 'l');
% Ec3
EcCom = changeRxnBounds(EcCom, {'Ec3IEX_arg__L[u]tr'; 'Ec3IEX_phe__L[u]tr'}, -exRate, 'l');
EcCom = changeRxnBounds(EcCom, {'Ec3IEX_met__L[u]tr'; 'Ec3IEX_lys__L[u]tr'}, 0, 'l');
% Ec4
EcCom = changeRxnBounds(EcCom, {'Ec4IEX_arg__L[u]tr'; 'Ec4IEX_phe__L[u]tr'}, 0, 'l');
EcCom = changeRxnBounds(EcCom, {'Ec4IEX_met__L[u]tr'; 'Ec4IEX_lys__L[u]tr'}, -exRate, 'l');
% allow production of anything for each member
EcCom.ub(EcCom.indCom.EXsp(:)) = 1000;
%% 
% Before the calculation, print the community uptake bounds for checking using 
% |printUptakeBoundCom|:

printUptakeBoundCom(EcCom, 1);
%% 
% Values under 'Comm.' are the community uptake bounds (+ve for uptake) and 
% values under 'Ec1' are the Ec1-specific uptake bounds (-ve for uptake). 
% 
% Create an option structure for calling SteadyCom and call the function. There 
% are a range of options available, including setting algorithmic parameters, 
% fixing growth rates for members, adding additional linear constraints in a general 
% format, e.g., for molecular crowding effect. See |help SteadyCom| for more options.

options = struct();
options.GRguess = 0.5;  % initial guess for max. growth rate
options.GRtol = 1e-6;  % tolerance for final growth rate
options.algorithm = 1;  % use the default algorithm (simple guessing for bounds, followed by matlab fzero)
[sol, result] = SteadyCom(EcCom, options);
%% 
% The algorithm is an iterative procedure to find the maximum biomass at a given 
% growth rate and to determine the maximum growth rate that is feasible for the 
% required total biomass (default 1 gdw). Here the algorithm used is the simple 
% guessing for find upper and lower bounds (Iter 1 to 4 in the output) followed 
% by Matlab |fzero| (starting from the line '|Func-count|') to locate the root. 
% The maximum growth rate calculated is 0.73599 /h, stored in |result.GRmax|. 
% 
% The biomass for each organism (in gdw) is given by |result.BM|:

for jSp = 1:4
    fprintf('X_%s:  %.6f\n', EcCom.infoCom.spAbbr{jSp}, result.BM(jSp));
end
disp(result);
%% 
% |result.vBM| contains the biomass production rates (in gdw / h), equal to 
% |result.BM * result.GRmax.| Since the total community biomass is defaulted to 
% be 1 gdw, the biomass for each organism coincides with its relative abundance. 
% Note that the community uptake bounds in this sense are normalized per gdw of 
% the community biomass. So the lower bound for the exchange reaction |EX_glc__D[u]| 
% being 8 can be interpreted as the maximum amount of glucose available to the 
% community being at a rate of 8 mmol per hour for 1 gdw of community biomass. 
% Similarly, all fluxes in |result.flux| ($V^k_j$) has the unit mmol / h / [gdw 
% of comm. biomass]. It differs from the specific rate (traditionally denoted 
% by $v^k_j$) of an organism in the usual sense (in the unit of mmol / h / [gdw 
% of organism biomass]) by $V^k_j=X^kv^k_j$ where $X^k$ is the biomass of the 
% organism. |result.Ut| and |result.Ex| are the community uptake and export rates 
% respectively, corresponding to the exchange reactions in |EcCom.infoCom.EXcom|. 
% 
% |result.iter0| is the info for solving the model at zero growth rate and |result.iter| 
% records the info during iteration of the algorithm:


iter = [0, result.iter0, NaN; result.iter];
for j = 0 : size(iter, 1)
    if j == 0
        fprintf('#iter\tgrowth rate (mu)\tmax. biomass (sum(X))\tmu * sum(X)\tmax. infeasibility\tguess method\n');
    else
        fprintf('%5d\t%16.6f\t%21.6f\t%11.6f\t%18.6e\t%d\n', iter(j,:))
    end
end
%% 
% |mu * sum(X)| in the forth column is equal to the biomass production rate. 
% 
% The fifth column contains the maximum infeasibility of the solutions in each 
% iteration.
% 
% Guess method in the last column represents the method used for guessing the 
% growth rate solved in the current iteration:
% 
% 0: the default simple guess by $\mu_{\textrm{next}} =\mu_{\textrm{current}} 
% \;\sum_{k=1}^K X_k^{\textrm{current}}$ (_K_ is the total number of organisms)
% 
% 1: bisection method
% 
% 2: bisection or at least 1% away from the bounds if the simple guess is too 
% close to the bounds (<1%)
% 
% 3. 1% away from the current growth rate if the simple guess is too close to 
% the current growth rate
% 
% From the table, we can see that at the growth rate 0.742726 (iter 4), the 
% max. biomass is 0, while at growth rate 0.735372, max. biomass = 1.0008 > 1. 
% Therefore we have both an lower and upper bound for the max. growth rate. Then 
% fzero is initiated to solve for the max. growth rate that gives max. biomass 
% >= 1.
% 
% Two other algorithms for the iterative procedure are also implemented: simple 
% guessing only and the bisection method. Compare their results with simple guessing 
% + matlab fzero run above:

options.algorithm = 2;  % use the simple guessing algorithm
[sol2, result2] = SteadyCom(EcCom, options);
options.algorithm = 3;  % use the bisection algorithm
[sol3, result3] = SteadyCom(EcCom, options);
%% 
% The time used for each algorithm in the tested machine is:
% 
% (1) simple guess for bounds followed by Matlab fzero: 18 sec
% 
% (2) simple guess alone: 35 sec
% 
% (3) bisection: 70 sec
% 
% Algorithm (1) appears to be the fastest in most case although the simple guess 
% algorithm can sometimes also outperform it. The most conservative bisection 
% method can already guarantee convergence within around 20 iterations, i.e., 
% solving ~20 LPs for an optimality gap (|options.GRtol|) of 1e-6.
% 
% 
%% Analyzing Flux Variability Using SteadyComFVA
% Now we want to analyze the variability of the organism abundance at various 
% growth rates. Choose more options and call |SteadyComFVA|:

% percentage of maximum total biomass of the community required. 100 for sum(biomass) = 1 (1 is the default total biomass)
options.optBMpercent = 100;  
n = size(EcCom.S, 2);  % number of reactions in the model
% options.rxnNameList is the list of reactions subject to FVA. Can be reaction names or indices.
% Use n + j for the biomass variable of the j-th organism. Alternatively, use {'X_j'} 
% for biomass variable of the j-th organism or {'X_Ec1'} for Ec1 (the abbreviation in EcCom.infoCom.spAbbr)
options.rxnNameList = {'X_Ec1'; 'X_Ec2'; 'X_Ec3'; 'X_Ec4'};
options.optGRpercent = [89:0.2:99, 99.1:0.1:100];  % perform FVA at various percentages of the maximum growth rate, 89, 89.1, 89.2, ..., 100
[fvaComMin,fvaComMax] = SteadyComFVA(EcCom, options);
%% 
% Similar to the output by |fluxVariability|, |fvaComMin| contains the minimum 
% fluxes corresponding to the reactions in |options.rxnNameList|. |fvaComMax| 
% contains the maximum fluxes. options.rxnNameList can be supplied as a (#rxns 
% + #organism)-by-K matrix to analyze the variability of the K linear combinations 
% of flux/biomass variables in the columns of the matrix. See |help SteadyComFVA| 
% for more details.
% 
% We would also like to compare the results against the direct use of FBA and 
% FVA by calling |optimizeCbModel| and |fluxVariability|:

optGRpercentFBA = [89:2:99 99.1:0.1:100];  % less dense interval to save time because the results are always the same for < 99%
nGr = numel(optGRpercentFBA);
[fvaFBAMin, fvaFBAMax] = deal(zeros(numel(options.rxnNameList), nGr));
% change the objective function to the sum of all biomass reactions
EcCom.c(:) = 0;
EcCom.c(EcCom.indCom.spBm) = 1;
EcCom.csense = char('E' * ones(1, numel(EcCom.mets)));
s = optimizeCbModel(EcCom);  % run FBA
grFBA = s.f;
for jGr = 1:nGr
    fprintf('Growth rate %.4f :\n', grFBA * optGRpercentFBA(jGr)/100);
    [fvaFBAMin(:, jGr), fvaFBAMax(:, jGr)] = fluxVariability(EcCom, optGRpercentFBA(jGr), 'max', EcCom.infoCom.spBm, 2);
end
%% 
% Plot the results to visualize the difference (see also Figure 2 in ref. [1]):

grComV = result.GRmax * options.optGRpercent / 100;  % vector of growth rates tested
lgLabel = {'{\itEc1 }';'{\itEc2 }';'{\itEc3 }';'{\itEc4 }'};
col = [235 135 255; 0 235 0; 255 0 0; 95 135 255 ]/255;  % color
f = figure;
% SteadyCom
subplot(2, 1, 1);
hold on
x = [grComV(:); flipud(grComV(:))];
for j = 1:4
    y = [fvaComMin(j, :), fliplr(fvaComMax(j, :))];
    p(j, 1) = plot(x(~isnan(y)), y(~isnan(y)), 'LineWidth', 2);
    p(j, 1).Color = col(j, :);
end
tl(1) = title('\underline{SteadyCom}', 'Interpreter', 'latex');
tl(1).Position = [0.7 1.01 0];
ax(1) = gca;
ax(1).XTick = 0.66:0.02:0.74;
ax(1).YTick = 0:0.2:1;
xlim([0.66 0.74])
ylim([0 1])

lg = legend(lgLabel);
lg.Box = 'off';
yl(1) = ylabel('Relative abundance');
xl(1) = xlabel('Community growth rate (h^{-1})');
% FBA
grFBAV = grFBA * optGRpercentFBA / 100;
x = [grFBAV(:); flipud(grFBAV(:))];
subplot(2, 1, 2);
hold on
% plot j=1:2 only because 3:4 overlap with 1:2
for j = 1:2
    y = [fvaFBAMin(j, :), fliplr(fvaFBAMax(j, :))] ./ x';
    % it is possible some values > 1 because the total biomass produced is
    % only bounded below when calling fluxVariability. Would be strictly
    % equal to 1 if sum(biomass) = optGRpercentFBA(jGr) * grFBA is constrained. Treat them as 1.
    y(y>1) = 1;
    p(j, 2)= plot(x(~isnan(y)), y(~isnan(y)), 'LineWidth', 2);
    p(j, 2).Color = col(j, :);
end
tl(2) = title('\underline{Joint FBA}', 'Interpreter', 'latex');
tl(2).Position = [0.55 1.01 0];
ax(2) = gca;
ax(2).XTick = 0.52:0.02:0.58;
ax(2).YTick = 0:0.2:1;
xlim([0.52 0.58])
ylim([0 1])
xl(2) = xlabel('Community growth rate (h^{-1})');
yl(2) = ylabel('Relative abundance');
ax(1).Position = [0.1 0.6 0.5 0.32];
ax(2).Position = [0.1 0.1 0.5 0.32];
lg.Position = [0.65 0.65 0.1 0.27];
%% 
% The direct use of FVA compared to FVA under the SteadyCom framework gives 
% very little information on the organism's abundance. The ranges for almost all 
% growth rates span from 0 to 1. In contrast, |SteadyComFVA| returns results with 
% the expected co-existence of all four mutants. When the growth rates get closer 
% to the maximum, the ranges shrink to unique values.
% 
% 
%% Analyze Pairwise Relationship Using SteadyComPOA
% Now we would like to see at a given growth rate, how the abundance of an organism 
% influences the abundance of another organism. We check this by iteratively fixing 
% the abundance of an organism at a level (independent variable) and optimizing 
% for the maximum and minimum allowable abundance of another organism (dependent 
% variable). This is what |SteadyComPOA| does.
% 
% Set up the option structure and call |SteadyComPOA|. |Nstep| is an important 
% parameter to designate how many intermediate steps are used or which values 
% between the min and max values of the independent variable are used for optimizing 
% the dependent variable. |savePOA| options must be supplied with a non-empty 
% string or a default name will be used for saving the POA results. By default, 
% the function analyzes all possible pairs in |options.rxnNameList|. To analyze 
% only particular pairs, use |options.pairList|. See |help SteadyComPOA| for more 
% details.

options.savePOA = ['POA' filesep 'EcCom'];  % directory and fila name for saving POA results
options.optGRpercent = [99 90 70 50];  % analyze at these percentages of max. growth rate
% Nstep is the number of intermediate steps that the independent variable will take different values
% or directly the vector of values, e.g. Nsetp = [0, 0.5, 1] implies fixing the independent variable at the minimum,
% 50% from the min to the max and the maximum value respectively to find the attainable range of the dependent variable.
% Here use small step sizes when getting close to either ends of the flux range
a = 0.001*(1000.^((0:14)/14));
options.Nstep = sort([a (1-a)]);
[POAtable, fluxRange, Stat, GRvector] = SteadyComPOA(EcCom, options);
%% 
% POAtable is a _n_-by-_n_ cell if there are _n_ targets in |options.rxnNameList|. 
% |POAtable{i, i}| is a _Nstep_-by-1-by-_Ngr_ matrix where _Nstep_ is the number 
% of intermediate steps detemined by |options.Nstep| and _Ngr_ is the number of 
% growth rates analyzed. |POAtable{i, i}(:, :, k)| is the values at which the 
% _i_-th target is fixed for the community growing at the growth rate |GRvector(k)|. 
% POAtable{i, j} is a _Nstep_-by-2-by-_Ngr_ matrix where |POAtable{i, j}(:, 1, 
% k)| and |POAtable{i, j}(:, 2, k)| are respectively the min. and max. values 
% of the _j_-th target when fixing the _i_-th target at the corresponding values 
% in |POAtable{i, i}(:, :, k)|. |fluxRange| contains the min. and max. values 
% for each target (found by calling |SteadyComFVA|). |Stat| is a _n_-by-_n-by-Ngr_ 
% structure array, each containing two fields: |*.cor|, the correlatiion coefficient 
% between the max/min values of the dependent variable and the independent variable, 
% and |*.r2|, the R-squred of linear regression. They are also outputed in the 
% command window during computation. All the computed results are also saved in 
% the folder 'POA' starting with the name 'EcCom', followed by 'GRxxxx' denoting 
% the growth rate at which the analysis is performed.
% 
% Plot the results (see also Figure 3 in ref. [1]):

nSp = 4;
spLab = {'{\it Ec1 }';'{\it Ec2 }';'{\it Ec3 }';'{\it Ec4 }'};
mark = {'A', 'B', 'D', 'C', 'E', 'F'};
nPlot = 0;
for j = 1:nSp
    for k = 1:nSp
        if k > j
            nPlot = nPlot + 1;
            ax(j, k) = subplot(nSp-1, nSp-1, (k - 2) * (nSp - 1) + j);
            hold on
            for p = 1:size(POAtable{1, 1}, 3)
                x = [POAtable{j, j}(:, :, p);POAtable{j, j}(end:-1:1, :, p);...
                    POAtable{j, j}(1, 1, p)];
                y = [POAtable{j, k}(:, 1, p);POAtable{j, k}(end:-1:1, 2, p);...
                        POAtable{j, k}(1, 1, p)];
                plot(x(~isnan(y)), y(~isnan(y)), 'LineWidth', 2)
            end
            xlim([0.001 1])
            ylim([0.001 1])
            ax(j, k).XScale = 'log';
            ax(j, k).YScale = 'log';
            ax(j, k).XTick = [0.001 0.01 0.1 1];
            ax(j, k).YTick = [0.001 0.01 0.1 1];
            ax(j, k).YAxis.MinorTickValues=[];
            ax(j, k).XAxis.MinorTickValues=[];
            ax(j, k).TickLength = [0.03 0.01];
            xlabel(spLab{j});
            ylabel(spLab{k});
            tx(j, k) = text(10^(-5), 10^(0.1), mark{nPlot}, 'FontSize', 12, 'FontWeight', 'bold');
        end
    end
end
lg = legend(strcat(strtrim(cellstr(num2str(options.optGRpercent(:)))), '%'));
lg.Position = [0.7246 0.6380 0.1700 0.2015];
lg.Box='off';
subplot(3, 3, 3, 'visible', 'off');
t = text(0.2, 0.8, {'% maximum';'growth rate'});
for j = 1:nSp
    for k = 1:nSp
        if k>j
            ax(j, k).Position = [0.15 + (j - 1) * 0.3, 0.8 - (k - 2) * 0.3, 0.16, 0.17];
            ax(j, k).Color = 'none';
        end
    end
end
%% 
% There are two patterns observed. The two pairs showing negative correlations, 
% namely Ec1 vs Ec4 (panel D) and Ec2 vs Ec3 (panel C) are indeed competing for 
% the same amino acids with each other (Ec1 and Ec4 competing for Lys and Met; 
% Ec2 and Ec4 competing for Arg and Phe). Each of the other pairs showing positive 
% correlations are indeed the cross feeding pairs, e.g., Ec1 and Ec2 (panel A) 
% cross feeding on Arg and Lys. See ref. [1] for more detailed discussion.
% 
% 
%% Parallelization and Timing
% SteadyCom in general can be finished within 20 iterations, i.e. solving 20 
% LPs (usually faster if using Matlab |fzero|) for an accuracy of 1e-6 for the 
% maximum community growth rate. The actual computation time depends on the size 
% of the community metabolic network. The current |EcCom| model has 6971 metabolites 
% and 9831 reactions. It took 18 seconds for a MacBook Pro with 2.5 GHz Intel 
% Core i5, 4 GB memory running Matlab R2016b and Cplex 12.7.1.
% 
% Since the FVA and POA analysis can be time-consuming for large models with 
% a large number of reactions to be analyzed, SteadyComFVA and SteadyComPOA support 
% parrallelization using the Matlab Distributed Computing Toolbox (|parfor| for 
% SteadyComFVA and |spmd| for SteadyComPOA). 
% 
% Test SteadyComFVA with 2 threads:

options.rxnNameList = EcCom.rxns(1:100);  % test FVA for the first 50 reactions
options.optGRpercent = 99;
options.algorithm = 1;
options.threads = 1;  % test single-thread computation first
options.verbFlag = 0;  % no verbose output
tic;
[minF1, maxF1] = SteadyComFVA(EcCom, options);  
t1 = toc;
if isempty(gcp('nocreate'))
    parpool(2);  % start a parallel pool
end
options.threads = 2;  % two threads (0 to use all available workers)
tic;
[minF2, maxF2] = SteadyComFVA(EcCom, options);  % test single-thread computation first
t2 = toc;
fprintf('Maximum difference between the two solutions: %.4e\n', max(max(abs(minF1 - minF2)), max(abs(maxF1 - maxF2))));
fprintf('\nSingle-thread computation: %.0f sec\nTwo-thread computation: %.0f sec\n', t1, t2);
%% 
% If there are many reactions to be analyzed, use |options.saveFVA| to give 
% a relative path for saving the intermediate results. Even though the computation 
% is interrupted, by calling |SteadyComFVA| with the same |options.saveFVA|, the 
% program will detect previously saved results and continued from there.
% 
% Test SteadyComPOA with 2 threads:

options.rxnNameList = EcCom.rxns(find(abs(result.flux) > 1e-2, 6));
options.savePOA = 'POA/EcComParallel';  % save with a new name
options.verbFlag = 3;
options.threads = 2;
options.Nstep = 5;  % use a smaller number of steps for test
tic;
[POAtable1, fluxRange1] = SteadyComPOA(EcCom, options);
t3 = toc;
%% 
% The parallelization code uses |spmd| and will redistribute jobs once any of 
% the workers has finished to maximize the computational efficiency.

options.savePOA = 'POA/EcComSingeThread';  
options.threads = 1;
tic;
[POAtable2, fluxRange2] = SteadyComPOA(EcCom, options);
t4 = toc;
dev = 0;
for i = 1:size(POAtable1, 1)
    for j = i:size(POAtable1, 2)
        dev = max(max(max(abs(POAtable1{i, j} - POAtable2{i, j}))));
        dev = max(dev, max(max(abs(fluxRange1 - fluxRange2))));
    end
end
fprintf('Maximum difference between the two solutions: %.4e\n', dev);
fprintf('\nSingle-thread computation: %.0f sec\nTwo-thread computation: %.0f sec\n', t4, t3);
%% 
% The advantage will be more significant for more targets to analyzed and more 
% threads used. Similar to |SteadyComFVA|, |SteadyComPOA| also supports continuation 
% from previously interrupted computation by calling with the same |options.savePOA|.
% 
% 
%% REFERENCES
% _[1]_ Chan SHJ, Simons MN, Maranas CD (2017) SteadyCom: Predicting microbial 
% abundances while ensuring community stability. PLoS Comput Biol 13(5): e1005539. 
% https://doi.org/10.1371/journal.pcbi.1005539
% 
% _[2]_ Khandelwal RA, Olivier BG, Röling WFM, Teusink B, Bruggeman FJ (2013) 
% Community Flux Balance Analysis for Microbial Consortia at Balanced Growth. 
% PLoS ONE 8(5): e64567. https://doi.org/10.1371/journal.pone.0064567
% 
%
##### SOURCE END #####
-->
</div></body></html>